diff --git a/single-nuclei-label.tif b/single-nuclei-label.tif new file mode 100644 index 0000000..0fee299 Binary files /dev/null and b/single-nuclei-label.tif differ diff --git a/slice-interpolation-for-zarpaint.ipynb b/slice-interpolation-for-zarpaint.ipynb new file mode 100644 index 0000000..51806ac --- /dev/null +++ b/slice-interpolation-for-zarpaint.ipynb @@ -0,0 +1,573 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5eed2e4f", + "metadata": {}, + "source": [ + "# Slice interpolation for labels in zarpaint" + ] + }, + { + "cell_type": "markdown", + "id": "856feb6e", + "metadata": {}, + "source": [ + "## Resources\n", + "\n", + "* Paper: Raya and Udupa \"Shape-Based Interpolation of Multidimensional Objects\" IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 9, NO. I , MARCH 1990\n", + "* [This StackOverflow post](https://stackoverflow.com/questions/48818373/interpolate-between-two-images) (not well executed for nd, but useful to read anyway)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "4048a74a", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import tifffile" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "a6d18acc", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.ndimage import distance_transform_edt\n", + "from scipy.interpolate import interpn" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "a2c32202", + "metadata": {}, + "outputs": [], + "source": [ + "import napari" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "1ac93484", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/genevieb/mambaforge/envs/napari-empanada/lib/python3.9/site-packages/napari/plugins/_plugin_manager.py:511: UserWarning: Plugin 'napari_skimage_regionprops2' has already registered a function widget 'duplicate current frame' which has now been overwritten\n", + " warn(message=warn_message)\n", + "/Users/genevieb/mambaforge/envs/napari-empanada/lib/python3.9/site-packages/napari_tools_menu/__init__.py:168: FutureWarning: Public access to Window.qt_viewer is deprecated and will be removed in\n", + "v0.5.0. It is considered an \"implementation detail\" of the napari\n", + "application, not part of the napari viewer model. If your use case\n", + "requires access to qt_viewer, please open an issue to discuss.\n", + " self.tools_menu = ToolsMenu(self, self.qt_viewer.viewer)\n" + ] + } + ], + "source": [ + "viewer = napari.Viewer()\n", + "viewer.open_sample('napari', 'cells3d')\n", + "membrane_layer = viewer.layers[0]\n", + "nuclei_layer = viewer.layers[1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24cb04eb", + "metadata": {}, + "outputs": [], + "source": [ + "single_nuclei_label = tifffile.imread('single-nuclei-label.tif')\n", + "viewer.add_labels(single_nuclei_label)\n", + "label_layer = viewer.layers[-1]\n", + "labels = label_layer.data" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "86393974", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAADGCAYAAAD7a3V2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAY20lEQVR4nO3de5SddX3v8fd39kwmd0i4hJALBAg0odWgMSAqRTkK5ViDVlmh1YWUY6DFC5WjgPWcenoOLVpBbXtA4pKClmVkCSxSpSqmqEURxIiYEAPhIhkSEpRLbmYyl+/5YzYnMzpJJjPPvsze79daWbPnt5/9PN/fzvOZ+c6zn/3syEwkSZI0Mi21LkCSJKkR2FRJkiQVwKZKkiSpADZVkiRJBbCpkiRJKoBNlSRJUgEq1lRFxFkRsS4i1kfEFZXajjRamAlpD/OgRhSVuE5VRJSAR4E3Ax3Aj4HzMvORwjcmjQJmQtrDPKhRVepI1SJgfWY+kZm7geXA4gptSxoNzIS0h3lQQ6pUUzUD2NDv+47ymNSszIS0h3lQQ2qt0HpjkLEBrzNGxFJgKUCJ0qvHM7lCpUhDt4sd7M7OwfbfkTITGpUqlIn95gEGZiLGjHl127TDCy5DOnDdzz9Pz/Ydg2aiUk1VBzCr3/czgY39F8jMZcAygMkxNU+OMypUijR09+fKSq3aTGhUqlAm9psHGJiJ9tmz8sjLLq1ELdIB2XjNZ/d6X6Ve/vsxMDci5kTEGGAJsKJC25JGAzMh7WEe1JAqcqQqM7sj4v3At4AScGNmrqnEtqTRwExIe5gHNapKvfxHZt4F3FWp9UujjZmQ9jAPakReUV2SJKkANlWSJEkFsKmSJEkqgE2VJElSAWyqJEmSCmBTJUmSVACbKkmSpALYVEmSJBXApkqSJKkANlWSJEkFsKmSJEkqgE2VJElSAWyqJEmSCmBTJUmSVACbKkmSpALYVEmSJBXApkqSJKkANlWSJEkFsKmSJEkqgE2VJElSAWyqJEmSCmBTJUmSVACbKkmSpAK0juTBEfEUsA3oAbozc2FETAW+ChwNPAWcm5kvjKxMaXQwE9JAZkLNpIgjVW/MzAWZubD8/RXAysycC6wsfy81EzMhDWQm1BRGdKRqLxYDp5dv3wx8F7i8AtsZ9VrGj2fj+xaw69Ttg96fCcd9fBu5+VcA9O7YCb091SxRxTATQ9VS4tkPnLzXTLzssFvHMenbj5iJ0ctMDFG2QPTuawGY8f1eXjiulZ0z9rWgqmGkTVUC346IBG7IzGXAtMzcBJCZmyLi8JEW2Yjyta/k7TfezbsmfocppfF7Xe7xldvZnX0HFP/4ax/m2I/+2F8i9c1MDNcpr+DJxRO4/92f3mcmADa9djsv9raYidHBTIzAUfM38fTq6ftc5tlTSvSWskoVaV9G2lS9LjM3lgNxd0T8YqgPjIilwFKAsez7B2gjaZkwgV/84zyuPPUulh60EfYz92PbJv7/2z9Z8hnOPelddP39NNq+85MKV6phMhMH6EAzATC9dSLT6cvEyTs/zLHXPUH3s5srXquGpZBMlKZMqVR9dW1/DRUBPe0DG6qxs7exc+tYWl5sq2BlGsyIzqnKzI3lr1uAO4BFwOaImA5Q/rplL49dlpkLM3NhG+0jKWPUaJk0iY23zGb9WcvKvzwOzEEt4/jWvK/zP2/4F7resnD/D1DVmYkDU0Qmfn7hP/PhH3zHTNSpojJRmjihWiWPers2TLKhqpFhN1URMSEiJr18G3gLsBpYAZxfXux84M6RFtkISlOm8PgX5rDqNbdQipG9P+D0cb189Lov0/VfXl1QdSqCmTgwRWWiLUqcMa7HTNQhM1EjvhJYMyP57T4NuDcifgY8AHwjM78JXA28OSIeA95c/r6pRXs7Hf9yBI+e9qURN1QvO2t8J1d8/kv+dV5fzMQQmYmmYSbUVIZ9TlVmPgG8cpDxXwNnjKSoRpMnncBdr7oOmLjfZQ/EW8Z3Meb6G/nbiy7wHKs6YCaGrmXcWO446QtUIhMXvTM44XvtZGdnoevWgTMTQ9c7uZuWba0VOcp0xLwtbH5+Mrl5bPEr1wBeUb0KOj7Sw8zWYn95vOz0cb188Pqv0v0mX/bQ6PHMe0/k0FKpIutec/b/pWX2jIqsW6qU8Qf/hmypzOt2m9YdzucWLadl2q6KrF972FQ1gHMmbOfJP2klWitx2TGpeN1veImDWsZVZN3jW8Yw+aYXaT3m6IqsX6qEXU9PInqiIuuOXrj84T/hhCM3Q2U2oTKbqgaxevE/8fSVi2pdhlQXls/5DyZ9eSst45vn0hTSvvzml5N44leHcMwfPFPrUhqaTVWFtc45itfPeqLi2xnfMobT37bKv86lsr+fuYJob45LU2h0653QQ+/BXRXfzq7NE3jvzB/SO8GL5VaKTVWFvbjwCG6YeV9VtnXdjB/RftMOaKnMuSrSaHJkazuPXf57tS5D2q+Jh+3gqBm/rvh2oiv4/kvH83/+8HZfBqwQm6oG8w9H3cHOc3xLudQebYw94aValyHVj4S1LxzB2eM3MHbWtlpX05BsqhrMsW0TeeHd22mZ4NWHpSMmb6N08EG1LkOqG9s7x9AWLfzZ3AfJVq8SWjSbqga0+pRb2H2KL3tId8/7N148c16ty5DqxkuPT+HOHTP42KHraD3USywUzaaqwiZ/8xFe//A7qr7d2Vc9WvVtSkN11Md383T39qps68W5/phTfdvRMYkFh3SQY6pz5OihHbMB+NTC26qyvWbiT5sK6922jRd2VOZ6PPvy34/4Nlv/9JSqb1cait71v+T0ez5IV1b+XUgrl36KLX95KoRn5qo+RXcwpXUn8098uionkN/2w0Vc+/wxnDp2MxPneN5hkWyqqqBn9UHs7N1d1W2eOGYcO899idLkyVXdrjQU2bWbEy56hPOeOLPi25reOpG7r/wHel+/oOLbkobrpy/OYvlxdzD+qK0V31Z0Bf98z5t5rGsc7z7uAbLNc6uKYlNVBcdcs5oHOqv/mUurXnMLOfvIqm9XGoreXbvo2HZwVbZ1aGkCvW3+uFP9evjnR/N8bzcT2qvzB3h0B13ZykemPk6O97pVRfGnTBX0bN3KpZ+7uNZlSHXnkPft5KKO11ZlWx1vHFOV7UjDEV3Bu9e+hxV/cDNjZu6oyjZvf8HPjC2aTVWVzLjjaS7d5PWjpP66n9nIhj89gvdteF3Ft/X+d9xV8W1II7HhycNY1TmVuxZdX5XG6htrf7/i22g2NlVV0r2hg7tvW1SVE3Ol0aRn/ZNsfOeUqh2xkupVy29a+McNZzC7dTz/cXJ1GisVy6aqimZ96gFOuOMva12GVHe6N3TQ8a5Dq3LESqpn6342m4s73tD3Bgsbq1HHpqqKsrubY77WxW3bfUee9Nu6f7mBjecdxkUdr/WIrppXwj2Pz2Xt7p3MbJ3It15urLwiyKhgU1Vlpe+u4oYL3sHybVNqXYpUd7qfeIoN/3Uib7j8En6wq7fQdd/yy9cUuj6pUnLzWN567yWs2f0bZrdO5Ien3MD/OOt2egt+l95hh/j5f0WzqaqB+MFDfHHpOazYMb7WpUh1p+e55zjoX3/ElX91McfdcwFbeop5+WPCp/0MQI0iz7Wz+Id/waNdO5hSGs97J2/hM2/6CnNfuaGwz+y7+oTbC1mP9rCpqpGW7/2Uf7rgXG7aenjFtvHJX8+jZZuvx2t0GnfnAxz354/w3jMv4H89N39E61r20pGM2WIWNLrk5rGcfe/7eXh332f0nTNhO/92wgq++bZriWkj+9y+nLqbaaXtfUeEu3xtsSg2VTXUcu9D3PLfzuYbO4u/MGhndvHl286g+5cbCl+3VC3Z2UnPI4/yo/e8go9tfgWbhvF5gbdtn8zX/vzN9K7+RQUqlCort7Tzjh/8BY939e37bVHi+LYJ3Hnq9ZSO+M2wrobee1A3y0+7gePbxvLJDWfTsrNUdNlNy6aqxlrufYjPXnAex3z7QpY8+abC1vt7376Y2f/7/sLWJ9VS78O/YNWisbztEx/hya6hNVbbe3dx3D0XcMP5byfu+1mFK5QqJ7e0c+a9H+Ajz57Erdv7XsY+ccw41px2Ix9909eH/EHM2QJzX7mBr77xeha1t/HRZxey5qdHV7Dy5tNa6wIELf/5U+b+J2ybOYMzb34r35r39WGvqzO7mHf3xcz70GP09PoOKjWO7NrN1BvvY+njH2D3pDbedNW9HDd2M0smPkcpBv59uHzbFP72X8/j2KvuB3OgBpBb2rl9y8l8bdxr2PaHK7jwoGdpixIXH/wM8//oBl7sGc+Hf3wuPZ0lWl5s+53H9x7cxTtfuYqrp/0EKPGRZ0/ijvteQ/ixf4WKzH0/oxFxI/BWYEtm/n55bCrwVeBo4Cng3Mx8oXzflcCFQA/wwcz81v6KmBxT8+Q4Y/izaCClQw/h0SvmcuM7Ps9pw3hV8NiVF3Dc+Q/7i2SY7s+VbM3n93mCgZmoH6XJk3ns4yfS2+/E3egJjv+7X9Dzwgs1rKxx1Esm2mfPyiMvu3T4E2kgve29vG7Bo/zdzK8zu3XigPu29+7ir545g+7c84dGC8m1M+/moJZxAFy+eQFf+94pYEM1LBuv+SydT28YNBNDaapOA7YDX+oXlk8Bz2fm1RFxBTAlMy+PiPnAV4BFwJHAd4DjM/d90Rl/gfyu0ry59I4fw6zrnuILs36w3+W7sofj//0i5n34UXq2Vv5TzhvVEH+BmAk1jXrJhE3V7+qd3E209vLZU5fztgk797t8T/Zy2bOLWHHfqwlPTh+2fTVV+335LzO/HxFH/9bwYuD08u2bge8Cl5fHl2dmJ/BkRKynLzj3DavyJtaz9jEAnnn3Mcw/d/9XYS91wQmf+wk9nZ2VLq3pmQlpIDNRGy1b+36Ff+i7f8bnj3tmv8t39ZZ4fPUMothLwKmf4Z5TNS0zNwFk5qaIePm6ADOAH/VbrqM8pmHqeewJZl31xJCW9UhuTZkJaSAzUSUt20use2j2kJb1+FRlFX2i+mD/X4P+ro+IpcBSgLF4EUw1LDMhDTSsTJSm+CkUqn/DvaTC5oiYDlD+uqU83gHM6rfcTGDjYCvIzGWZuTAzF7bRPswypLphJqSBCs1EaeKEihYrFWG4TdUK4Pzy7fOBO/uNL4mI9oiYA8wFHhhZidKoYCakgcyEms5+X/6LiK/Qd7LhoRHRAfwNcDVwa0RcCDwNvAsgM9dExK3AI0A3cMn+3tEhjTZmQhrITEh9hvLuv/P2cteg7/fOzKuAq0ZSlFTPzIQ0kJmQ+vgxNZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFWC/TVVE3BgRWyJidb+xT0TEMxHxUPnf2f3uuzIi1kfEuog4s1KFS7ViJqSBzITUZyhHqm4Czhpk/DOZuaD87y6AiJgPLAFOLD/muogoFVWsVCduwkxI/d2EmZD231Rl5veB54e4vsXA8szszMwngfXAohHUJ9UdMyENZCakPiM5p+r9EfFw+bDvlPLYDGBDv2U6ymNSMzAT0kBmQk1luE3V9cCxwAJgE3BNeTwGWTYHW0FELI2IByPiwS46h1mGVDfMhDRQoZno2b6jIkVKRRpWU5WZmzOzJzN7gS+w59BtBzCr36IzgY17WceyzFyYmQvbaB9OGVLdMBPSQEVnojRxQmULlgowrKYqIqb3+/btwMvv+FgBLImI9oiYA8wFHhhZiVL9MxPSQGZCzah1fwtExFeA04FDI6ID+Bvg9IhYQN8h26eAiwAyc01E3Ao8AnQDl2RmT0Uql2rETEgDmQmpT2QO+lJ2VU2OqXlynFHrMiTuz5VszecHO+ejqsyE6kW9ZKJ99qw88rJLa12GxMZrPkvn0xsGzYRXVJckSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQD7baoiYlZE3BMRayNiTUR8qDw+NSLujojHyl+n9HvMlRGxPiLWRcSZlZyAVG1mQhrITEh9hnKkqhu4LDPnAacAl0TEfOAKYGVmzgVWlr+nfN8S4ETgLOC6iChVonipRsyENJCZkBhCU5WZmzJzVfn2NmAtMANYDNxcXuxm4Jzy7cXA8szszMwngfXAooLrlmrGTEgDmQmpzwGdUxURRwMnAfcD0zJzE/QFCji8vNgMYEO/h3WUx357XUsj4sGIeLCLzmGULtWemZAGqlQmerbvqGjdUhGG3FRFxETgNuDSzNy6r0UHGcvfGchclpkLM3NhG+1DLUOqG2ZCGqiSmShNnFBUmVLFDKmpiog2+oJyS2beXh7eHBHTy/dPB7aUxzuAWf0ePhPYWEy5Un0wE9JAZkIa2rv/AvgisDYzr+131wrg/PLt84E7+40viYj2iJgDzAUeKK5kqbbMhDSQmZD6tA5hmdcB7wF+HhEPlcc+BlwN3BoRFwJPA+8CyMw1EXEr8Ah97wi5JDN7ii5cqiEzIQ1kJiSG0FRl5r0M/vo3wBl7ecxVwFUjqEuqW2ZCGshMSH28orokSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQA2VZIkSQWwqZIkSSqATZUkSVIBbKokSZIKYFMlSZJUAJsqSZKkAthUSZIkFcCmSpIkqQD7baoiYlZE3BMRayNiTUR8qDz+iYh4JiIeKv87u99jroyI9RGxLiLOrOQEpGozE9JAZkLq0zqEZbqByzJzVURMAn4SEXeX7/tMZn66/8IRMR9YApwIHAl8JyKOz8yeIguXashMSAOZCYkhHKnKzE2Zuap8exuwFpixj4csBpZnZmdmPgmsBxYVUaxUD8yENJCZkPoc0DlVEXE0cBJwf3no/RHxcETcGBFTymMzgA39HtbBvsMljVpmQhrITKiZDbmpioiJwG3ApZm5FbgeOBZYAGwCrnl50UEenoOsb2lEPBgRD3bReaB1SzVnJqSBKpmJnu07KlO0VKAhNVUR0UZfUG7JzNsBMnNzZvZkZi/wBfYcuu0AZvV7+Exg42+vMzOXZebCzFzYRvtI5iBVnZmQBqp0JkoTJ1R2AlIBhvLuvwC+CKzNzGv7jU/vt9jbgdXl2yuAJRHRHhFzgLnAA8WVLNWWmZAGMhNSn6G8++91wHuAn0fEQ+WxjwHnRcQC+g7ZPgVcBJCZayLiVuAR+t4Rconv6FCDMRPSQGZCAiLzd17Grn4REc8BO4Bf1bqWGjmU5p071Nf8j8rMw2pdhJmoq32iFupp/maiPtTTPlEL9TT/vWaiLpoqgIh4MDMX1rqOWmjmuYPz35tmfl6aee7g/PemmZ+XZp47jJ75+zE1kiRJBbCpkiRJKkA9NVXLal1ADTXz3MH5700zPy/NPHdw/nvTzM9LM88dRsn86+acKkmSpNGsno5USZIkjVo1b6oi4qyIWBcR6yPiilrXUwnlz7zaEhGr+41NjYi7I+Kx8tcp/e67svx8rIuIM2tTdTEiYlZE3BMRayNiTUR8qDzeFPMfDjPR2PuEmThwZqKx94mGykRm1uwfUAIeB44BxgA/A+bXsqYKzfM04FXA6n5jnwKuKN++Avhk+fb88vPQDswpPz+lWs9hBHOfDryqfHsS8Gh5jk0x/2E8X2aiwfcJM3HAz5eZaPB9opEyUesjVYuA9Zn5RGbuBpYDi2tcU+Ey8/vA8781vBi4uXz7ZuCcfuPLM7MzM58E1rPn87JGnczclJmryre3AWvp+zT6ppj/MJiJPg27T5iJA2Ym+jTsPtFImah1UzUD2NDv+47yWDOYlpmboG+HAg4vjzfscxIRRwMnAffThPMfomaef9PtE2ZiSJp5/k23T4z2TNS6qYpBxpr97YgN+ZxExET6PsH+0szcuq9FBxkb9fM/AM0+/8E05HNiJoas2ec/mIZ8ThohE7VuqjqAWf2+nwlsrFEt1bb55U9wL3/dUh5vuOckItroC8otmXl7ebhp5n+Amnn+TbNPmIkD0szzb5p9olEyUeum6sfA3IiYExFjgCXAihrXVC0rgPPLt88H7uw3viQi2iNiDjAXeKAG9RUiIgL4IrA2M6/td1dTzH8YzESfht0nzMQBMxN9GnafaKhM1PpMeeBs+s70fxz461rXU6E5fgXYBHTR12FfCBwCrAQeK3+d2m/5vy4/H+uAP6p1/SOc++vpOyz7MPBQ+d/ZzTL/YT5nZqKB9wkzMaznzEw08D7RSJnwiuqSJEkFqPXLf5IkSQ3BpkqSJKkANlWSJEkFsKmSJEkqgE2VJElSAWyqJEmSCmBTJUmSVACbKkmSpAL8P8MHBpuBOIjlAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "first_slice = 30\n", + "last_slice = 38\n", + "image_1 = labels[first_slice, ...]\n", + "image_2 = labels[last_slice, ...]\n", + "\n", + "diff = labels[first_slice, ...] - labels[last_slice, ...]\n", + "\n", + "fig, ax = plt.subplots(ncols=3, figsize=(10,3))\n", + "ax[0].imshow(image_1)\n", + "ax[1].imshow(image_2)\n", + "ax[2].imshow(diff)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "b3c154f1", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.ndimage import distance_transform_edt\n", + "\n", + "def distance_transform(image):\n", + " \"\"\"Distance transform for a boolean image.\n", + " \n", + " Returns positive values inside the object,\n", + " and negative values outside.\n", + " \"\"\"\n", + " image = image.astype(bool)\n", + " edt = distance_transform_edt(image) - distance_transform_edt(~image)\n", + " return edt" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "b8271a7c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-202.5660386145713 -84.07044179766024 20.518284528683193\n", + "-206.0218435020908 -86.68854786664434 18.788294228055936\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAADGCAYAAAD7a3V2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACWBklEQVR4nO29W6xtTXYW9o2aa629zzn/+W99o/33D21QJ7GdgLnERCGKnKAEB0VxeCBqRwI/WOpIMQpEPGCTh+TFUhKByUOUKB0Z2ZEAxxIgLERCjAVCSAZjiBPcbmw3dsdu9+9ut7v93845e681a+ShalSNus3L2mufs87ZNXT2WXPWrHmfo75vXKqKmBldunTp0qVLly5dbibmWV9Aly5dunTp0qXLiyCdVHXp0qVLly5dupxAOqnq0qVLly5dunQ5gXRS1aVLly5dunTpcgLppKpLly5dunTp0uUE0klVly5dunTp0qXLCeTWSBURfRsR/RwRfY6Ivue2ztOly/MiXSe6dInS9aHLiyh0G+NUEdEA4OcB/HsAvgDgHwP4Dmb+2ZOfrEuX50C6TnTpEqXrQ5cXVW7LU/UtAD7HzL/IzNcAfhjAt9/Subp0eR6k60SXLlG6PnR5IeW2SNUbAH5FrX/Bl3Xpclel60SXLlG6PnR5IWVzS8elSlkSZySiTwH4FADQdvd7dx/6cLEn145UXeeZ7Y26AKioz9WLJ7Vv2K4qUn7c1vVU6i6V+pUhuaXakcN+XKtDyCPADABcnquIFOd1iu3peSa3c6g1uR0AaPI8bp1q2xYc4/q9r2J/9X7jQd9Iuk60rqdSd6l0nfB1nj+dmNUHINUJM+x+772HH8p25ca6HIoqdfLT1epM7ZOXqzLZpC8B2fLJr9HrIGe7cLonKO7uPtO1z3Hi/mqXvESytohr5fm2xvbJ8+dtaO0YU8fP5Prtr2J8VNeJ2yJVXwDwplr/GIAv6grM/GkAnwaAyzfe5Df/8//SPTgD9+s/ALfOodyV+RbCAEwcPxZZNrGMDIdyIgb54xFknd2vbCeGkTIAhjiW+zrGf72yLPXycvhyEaPeNhUt2LRw1lhb9dat3ya/rNZDmV9nX+Z+Xd1QhrjOfpu1cR1hm6wDbMmXu+2w/tdv18tky2WysZ4sU6gT6xIDsEh/5dXavH78jfXZrTfqkHXbf+b/+B9WvZcV0nWi60TXiSiz+gCkOvHSa2/y7/p3/ySIGVww/3VCng23jjO3vX1gp592IH+c+BzX2g1LrzGsj+5csP6cI8OMro4dCLwx4AGwG1Lkavn98QCwIdfuhG9SETr1nc7eQ9KeuWOy8aphKDxH1u1f0Q5SWC8IpTqPO3+6b3IcRTbzdabKcQB8/ge+v/mcbotU/WMAnyCirwfwqwA+CeA/nd2Lan+crKcP5bTgEYBjBjzmgCOAB2I5kILJalH7WiYM/tjMFI5rKQWOgRjGL1tfzwIYiMFMIN/wW2YQE0Zr1HkIAMMYpzfWEjjxCzjAIcNgCxDcdhgCLCN+hXExsRa8Xsuv08hS8eTeGQTy+BN+OX4PybG0gvj6mKrLXulv8n7mpesEuk64C4yLd1gnjtMHuOsSYlUjWDUg1/WmSJnelwJhWEHihJQSO7KgdFM/ylNfY1mnfHVM8k5RnKf1HKUsITjZcacIVfUaSZMlRaiMuy5NopjgDUVdRknbl583P0+LTMkzaZbnx0GsOyW3QqqY+UBEfwLA3wYwAPiLzPyZyZ3UhbNeL35TQKmBzFrwMMYmVnkABmMngQNAUq6BQ4NGYpnfsJESIJBlqOUBHADFkgcOIAES317CwiXUuV+CBWMw1lnivjy0rgCMSUHE4cICEDEMWK/Ihp01DkRQMPD7KlDxWJKCjAIY/YH702hgSUBCtpO/2glQmVOWm0jXia4TXSeiHKUPXhKwX0h2dL2pfZaSmkki54mVkAAIOY2fzkmvsfBYFWSK3PeXEYLFz5GE9MQ6jshwk1DJMYtnRZFQJR4pv171TGlSRRRJFurnzQnTpHdKtbHJcr6e152Q2/JUgZn/FoC/tXiHHBQA18jk4BFcgJw2JARlkU+ARwIWdfAYjF1kiWvg0BZ4AiaZhe7K7BFP1IlVfQu0Fa5BRcDEgKaBhNhb5gRYAytg4Y/PxLA2miYaRGLpShAhBoUWO7e21bqACBrAoH+RAYAGieQYrnGLoKL2E8v8FqXrRNeJrhNRVuuD7DdDfNZec+04UwQnJwt1rxMA6+rUvFWnvsZYmC0zCk9Pvv/cc9QEKBKoaUJVvUZKCVVCmiQEKNeof0mdX77R1jlJ3at67keRqXx7foyG3BqpWiX6prN12R6YPnLw4HS7B4YqeBBPgscgIFIBj0EBQ8sKz4FDQEMA46bWuQOKMS5LuV+y8BY5K9CoAIlJjufDH8aG0MdovcXOBGNsE0Ri490AEYo1QlyCPGwY8qZ/bm2rdRH1rmdDHqyUrgISQTEQLy+xzM9Fuk4skq4Td0gnvCwNx+VhMu3VqZGmNWG+fP8ixKUIh5AGYo7P/Lav0fqThO+E0Ar9Td7nEkKVSfW4VBKqqreKEEhgJFW03DulCGSVTKFSjrK8WA/LnNapyHmQKqACHlzeOCEm5GagkYCHD2OEMgiYtMMbOXhEy3zaEk/WFaAYsqVlnn2Bpvgy2iLgAJRWuKExAIplk4Y6BDQYSTkxwTDhYE0AEQl9DGY5iExa52KN52ACDyLB9ex1Nv8VyzwpzwEHBTjo70mDxBLL/KxQpOvEpHSduIM6gUgKZglBhTxMEYqpnKKpspa4ECBHMmPcA/U2wMmusRCVC+b2ychCdqzm/VAkO+54qBOqnE+uIFSpt0oty/NS7V2VTKnj10nQRFll/2Jdt72msm9FzopUhd/Q+FfWi7rxjwKIcFg/JXgMxhbAAbW9Bhx5qOPY/BHZPwGPDFTcr7PKUQMSJhjVSFoAG2ND6GPEcSAS214KJMA11ZSAvxN1AWHbDUIenP7m25ogUQGYcwOPrhPT0nXiDuqEkinPjpY8GXyqrt5eSyJvJZY3jydJ60MMAwYufaJrrBGrRKWS3nTlMavPUepLHlVLRedUl6YJVU6mkutc4J16qmRKhyMn9OJsSFVikaP1YDhpdFJw4QRoNEAIqBwDHhu13rLENXDkoY64HnNGhiNBBABG9TatNyFiyEKHOlIgOcCEbSQ+ViAJfbjnzj6XZL117nIzCTDsupVD3hfFHBJ52f596lySNSGPgAoaoDj7rYDErGV+RtJ1Ypl0ncCd0QktNUIz14NuiXerRlJqSdctL1ZSj1EkrbvOCX7bCa9RQo55Xd1e1K6x9hzDsBDyyda8VGsI1VAhVPJMVOjvKO+UydaT9lH9VvYt1sM1c3nsEJ7kSb04G1IFIAMEKcstdF0WfykAhyxH8HDdwuO2m4LHRlnnALAxzp+bAkoEDgGMm1rmuUU+eD+ygIpl43s/oQCSjd9HHS3klRysSXJKiPg465xJgYK2tH1LTerXX+fqkEfN+sbEtrWW+blJ14lJ6TpxB3ViQo7JjxKpeaOWlNfKEnInYUA1xlN85ae5xrQOIsESIlchns1nRIrc3JRQDalHquWhOoV3Khm3T7WX2jN8NJkyrAjgc0CqZi1yANWeTQE0oADDL0PKkIDGseAxKKu7Zolr4BALfFAgo8tl2xoZmSCeWMumCH0IoIxMCZBEse5LoaQEhmO4Azba52usc5KW2VjAGv+ufGquAchmSbrB/y3aI++1HfLQjX5eN9l2U8v8TKTrxLx0nVDvPKv7IupEcq+VT6VFFmrhrVadmrR6+dWSy1vXIYOstpLWb3qNLZFQFQ/LnhUIsAOm86i4PE568JJA2WGCTKnwZDxn5WYC8UFBnia9U9l2eS6hvuHME6WOkZOpwbcDE6zyLEgVgNDw540EsoexKMQBAYwUOPSo0EvBY0OZZa7AY2McIGxUTyaxwmtgApSJuEu6kkuibVj3yzmoBDAhjkAC1cvJW+gWJpQdrIk5JR4A3PEIg1+fAxFnmfuXmFjccCBCQDWXJDTeEyEPqaMtc6BqcZ/CMj8r6TrRlK4Td1MnWD0GV6C3tUN3+fY1ZEV7h/LjT+VUJWUM0JglrQtZOcE1JmI5XqchxCFXZq4RiuT4a24RqqbHjFJCZdUo7IFQac/VMd4pk61TpQ7K7e5+9f51MuWukdX1st+WtcUNORtStcoil1+aDnGkIIIAHuZI8NiQnbTEW8Ax1Y18iRjVXUR3G5cARW6RR4PXZ4coA9gfBWKhb4ydBJGpsAdJyxzCFQIWcso0STe03rJfAjZZyCNY0GWCbmz8T2+ZzynM05SuE23pOnE3dSIBR6UCAsRTidirT7Ui7LfkPOSnj5GkdVLv5ibXmK77BfX+FnvkKnlUNULVvphIqMTbxYb89DbpnxBLPW5W1ZaSe9DeKUWqlpCpcjlNPmd13ECmBi5DfQtfz3mQqgQUZHmBRe7L5kMcKXjoUaFPAR5bnz8SrXMNJrHn04DjAETETbPhlkcY6PF5YrfxaJG7/JC6hX6wWAQigAfwinUOIIwyLWBCHgwEdKq5JAtCHollri11IIKUwqKTWebnIl0nFknXCbmeO6ATQPDwAPE+5JprwF/zLM0lnE+Vt3K2lobqzOj6JtgBboBZ5a069hpbxIqDwtePmawTwJXE9BahKo5HkZQ5UhJJU/gbhCDROu+UmVmncp/a8uJ8qSPJlMh5kCoAwRWH1DpaYpED4duphjgSoPDWOFHsDr4WPDZhubTEW8Che0BJ2RoZEcMdzjKP4DEQMFLs2eSeSW6Rp+sONLAIRHLrXPQoZpCo9t63zknKQGj4XEijGvLw75rESq8AArFq6G/LMl/1Vm5Zuk5MSteJu6cTEjYFPBlF9sgoBek13p+lc+/lx23Nv1c9JshNdpxNTMwLkbt5jd7bBX0tFMnM1DUCkfi448KRPRuf5WQvRaoTqpCblXipZrxTcju5B0mTIE2aoJapvpyQKa/umkzxgDLEp64lewHtbV7OiFSpX33hoXzeItdlNYs8B5PbBg8BDlkGUot8WDE1h1FKIT2WBFQEUCwoAImBqVroezscBSLufN4a95a5e9y+AfYtc8glEbSRF1nrUq6RITAB9b4FRQh+PwrXfhuW+VmhB9B1Yka6TtwxnRCgDOsU78kzViGZgVwBxT3kpCTPmWr1rmuNH7U0PBjOOzq9sxsCDIGZkwFBb3KN+n5TQ6x9jSHsJ/tmhKq2fzyQIlRDTqxSMrXUO9UK9SXkStWvkqlG8nktxCf19XWU9zmzXcn5kCqg8sAqN+LBYsoih1qu5YzIclrOxfIceMj2pcAhgLFm1GiRENZADHcIqFgQRgk+sIWbh0T3bIoW+da4PJODHW4EIkxuwlkBMAl1hORZb4Xr1pyIwUxqu3qvoXGfscz9vqe0zGV5obH4dKXrRFO6TtwtnWAoUsixDAAS75U8m7ANaWVZrXibWkMNLE1mb4XtitCbJy1s4N6/jRfcGlR07hpdHVVvoJh03rpGQhr2s+yPgeR6as+sSqiETCWkKo6cXg0n+kNyjQBRPM9qMqUJ2TAT4lPXUTzQqe0VORtSlShwsRwt7jUWOYBqiKPouSSNfQEWNwePcFzcDEAADx5+eWSXPyKAYpi9VW7Ucmqha4scZlwMIpY45IswJNzBMHBj8EyFPAIA+cY7Wth+RUzlwgpHpYzjCZABjN+2yjIPy6qL8xlJ14l56Tpxt3QiAHC6Ok2uBMyRVZbVGc+TLpvzQi3NzyILmANj3Ma8I/FW1fKellyj+xBZEWYZZDO9Nr2vHdyfkDxihIFJ9bmaQy8MFH6ltx9M7PUn5Ebuec47VfxS3L9FplzdSvL5TUN8+fbnjVQljUbthlXDstYij2VpiCOxwpEDynLw2JhxMXCE5RVhDhEZLVrAA/AWORv/xRpEqxzILfStGRMQsRQVcApEwjLiGEBMnOSS5CGP0PNJrHH4Bk9Z0CFBN7OedXfyFCQUMMi51b7ayg7fS8Uyn1w+J+k6MStdJ+6WToTwkV7POKg8iLjO4Tm1yNVS79TSXn9LykjGrvIERHuraudrHbN2PdaY6Olp3YsnQ5FMcZ34ZBITzykjVggEZol3qgjrZZ6qpI58j5TX0UMfKBJWC/GdgkzRvEKcD6mCenCoPIDEMudq+ZxFLmU18EhDHOMi8JB6c+ARfqVLuf+61oCIZRNGi5b9LBs4w8BPxwE3r1luoe8xoAYic2EP40EiTDzrn5tO0m2FPPKeTzKyNMgPgEgOUDi8Q9UqJg17HsZIQWdVHsnkMmGJwjxt6TrRlq4TuHM6IdcYQkGYJ1jhmeqeAqoe5F7DOThZbnmfanlO7euuhN7Ye6t8Urdbb9/7kuR4YrjPeosk9FcL+9mN92L6oR7Cc2ncI5ARqo0K+XlCZQda553KQn3Vnn05+SKUyedCyvRQCPLrz5n8Fg+tsr1FptqvGcA5kar8ZsJDSIEkt7xby1MWOYDEAk9CH4V17sBjWwBJCh5bWW8AhwaNtQCiwWNEHJPHgjEGxPWWecVC32JsgsiUdS7WuDuau2bJ6fG6uzzkocCB5R1rQMgsdJ1HkoCFsr5X55EsXT4X6TrRlK4TUnbHdAIpvs0RLFLLhfeK43MXcuX2LW94qmfgmvn7igR0Hwa0W/LeHS56xM0mvgvB0vfgyU71GA1CNXdeIVR243OoNspDdYR3KniYJshUmldVyZeaCvH5cya/6Q2W21pEaqpeJudBqmpgUWxD+hBIbQ5WOFZZ5ABSwEAtxOFDGjcED/0LrOs+LuGREXo+MgcqBmIV2zBVh1joewyLQCS3zk0wd2KXdfZWnuSSrA55KCyAhD8IyopU73hpHsmcZc6qTLb7c7WWz0a6TkxK14k7qBNAfCZe5gjWEu9VUp4TgXBsKpa192eJR6t2PGKGGZW3x5bkpjYWVTNJ3XpSIVO/VK4x5FH5+lQ5Z3HPKtwXSRVCgvoi71RChFJCdRSZGhohPn++5Fc/JFS2tchUq+2dkPMgVUiVIQlzANlNxsZfL1O+LNVzqzv7S7fZAkyGbHkOPLY0NoEj6fV0zJg8yTg81vVuKqxxJMtLQURb52KNA26CWcADB98w5CGv1bf9hWXOusx9B5N5JHLMKQDieL4cXMrlBRrzFKXrxLR0nbh7OgGgCtjABMHKCFVkslIRyEODLXLlzuNJNZUJ5Ktvxb9jc2CMO4LdUMxtyutNiCMyHAmcids0CRRPkyNUaBKq5Nw63LchRaYQPGLB65U/N/8eYs+7jFhpMqXXVfJ5M8QnvzUS1CJTpyBS1PgwvJwNqQIQlT6sc/KQ9HdFpDapbaVljrCsQQWIFrksA0jCHUnSrS4/Ejxk35iku8IyV1oiXcWdBW5dTyZwsNABHAUiejk8F9/YBOA4IuQREnNDy+9auLQMsSXXIFADhgrotHJFEms7OV5l+Ryl60RTuk7cLZ0IjwOpwZGDOFQ9qTtLsDgb82qCXLVCcccMxSB1MBJoRCAuwzUXdVpjV8U6AEYGJJSYn5cAu/XH8GG/NYRq3KYeKvFcaYKmzxX2156pBWSqyJeSEJ8MJjr47TclUyckUlrOh1QtvcEKyNQsc5Mtt/JGYpnN6tk6kCwAjy2NTeAQ0BhWWOYjDAZvkbsJYx1suEZB8kQQQGPrr39vN4tAxIR8E4TlZG62I0MeDA3oPq9DAQH79ymJu808kpplngNEAgrz4AJ1nLA83fY9fek60ZSuE3dUJ7zkGBdIU459lNadJFisK2GWXOVyDKHS+w17Rpgrb0AyIOgUkQvix5gas3yqsHlDYegGslwMOFoc34f27MaRMUf4kHqnPKGa9U7lIT9NpoIXihPSJTlSjkzJtuwlLCVTVNlWlDeOIasLdeF8SBUwAxycPQwV8gDSZbXeBBWkeSOyXYOGKxMgiWAgfy1rfAo8BDjyXk8tcUm4qmcTSTm748H4MIiz0MEeNNjAZiGQFohsAez9Uxm8nW18ay25JEtCHmCCIWcsGeJwZ7llLm0YkrKJPJJamSoPnnsFKE1w8R9INdxxjtJ1opCuE42yF10nNLHJ9CInTsk+odI0wYqkSv4olLHPvarmC6FOoGrJ7K1R2eW8Ogw4VIZYmNpXQoB2a8rN3sNEFjBjnVAlg44Kodp6D9UmetF07lTTO9WaRDknUz5figcocgWXL1VLPM/OVSVLNaJ0KiK1QCnOi1TBP7wCONQqcXKjrTAHAMhI0YA7hLbM11jksR4Hi3yJNR7rpcAhoLFkSg49sKHu7VRa5PDAcAggssY6H0DQwCG5JGtCHkst86PySMSSJ4UjGlCg60+AwwSoFAp2JtJ1IpWuE10nqrw7POOsqn4Wqu4SgsVit3iC5Z47IokJ+5M/d71n4JJhFwDAHBypsFuARkwOsSDnFaJmRgsQuRHStRAw7tz1U4VQJcMmeFJjN4Rx60iV3cZcKpAiU9nzFHKUk6qETIXk9CxfSof4BgaGCpHy5znKK7WQSFVfT43QTbzG8yFVKxljUT0BDi7Wk1AH2sshtCHAogABgLKwnVXessanwCOWLwx1hHAE/DWYzCqXFgAFiIANYA5NEBn8r1EhjLUhDzAFT4gk6B5jmUeL27fwVSs8BZkAKPKNhBNybCUzIKmGOOYNkKcvXSea0nUCd1InmoQJKK93imRVcLJ6XKkrz5P9tEOWwzYhWW6/klDN35OQcUeQzAHggTDuyHmV7IK8LMugvQVvTRH6O1w4j+dwAEyDUMm9Su8+uyWMu0io2I+j5UKH+gBIQn16RHUhaPIb8qWyEN9k4rk6z2IypfevEa2sziJv1Arj4kakiog+D+BduOmvDsz8+4jodQD/O4CPA/g8gP+Emb+27ICV9ZxhZiw0D3ckmyvbRWo5JW7ZhkZzyCx2XW/wVriAQlwuwWNLhwI4wrg8C3JIBnXd0SIHYi8n9VwaINIKe1gi9RtzSZaEPKByR+TUoweElmWu5zkj/05lzrMkjyTco78t8l3RBVNqgIJ0WwEOFHepgcoprPKuE10nuk6kcnKd0Meu4F8zvwoIvLSoz+06yXMQAkUUvFjMXO35tnbsKlkmC5g943BpMO6A4apyG5XR2c3e4vr+xr9H/0365HKzd6FFub5i+hshVFtH5hyhcp4q77B1hCx7ToEs6RHVKS7HOkKcZDunhCrPlVLnqBKiJUSqWF5BpCqvipI6bYujDL6ul3+Hmb+ZmX+fX/8eAD/OzJ8A8ON+fbnUvrvihtN1mltXy6XFbhMgyfNGAFTyQFhtOw48DCy2dMCOxsk/B1QHbOng9kFaVvMExMThuF2HZ6I3gbM/G0ETyjuB6VybkKvjn72rmz7v+Offj+znFSQshwMF7U8VSL1MzpVH1eO8XlPRJrYdL10nuk6E59x1AsCpdIK5/qcvl8u/uH/2V6mf19HbxCsjg16GaVpkIExPJOK1tMF3SswBGK7ZeZmG+WOZawsww16YhCgdLlyvwmHPbR5APty3IxwuCYdLYLwk2J3bRgcfhtRE3IhXCzFEqNZ5A/CGYbfyJ/UY9sL98Y7BG64TKoJjJ5o0EbscK/k+SZUnZbXtomeeUNU+EL2fnFbp6VK5jfDftwP4Vr/8QwD+HoA/s/ooM0yyBA1kwBBzH9w6l+vqK0uWqcwb0fvoEIceuDAm6tbBY0djsMJ3PrC9pemguYEJ85sNNAar3EIsbG19q+fGQC2fxG3zIQyysExYGvJwvZpubpkDM5a53EdmUSdhjEq4Y7IesnpFXfXsTi9dJ7pOdJ1I5SidYEMZSeL0N5eQ61Q5lr5XvUter7FNiCn747v5+/wj9sMVcHa0pUMsMLlhFezWEaPtIy48Wlouf+krGD/4Muw25m4dLlyd4ToNISbkLBAqd57x0nmnxIGbkykgeqbK+f4Qk8/zRHVJPB84DJlQPFw5R0LouSzXZXq/fB+E11+UF/uE+vPkiWdcuDclVQzg/yJ3Jf8LM38awEeY+S13cn6LiD686EiESob/PGjkYrI6xTrK7bUwh6zneSNyjDxnRKz2OfAIVjlZ7DAqACpf5gjCwAzJLHTgscFAI9wwJutBJDT+3twwxJgLebhBetVz8vOoGeIb55GEVj0Agle4vBxohzuAAEJh1xo4KCCBApMT93TqOtF1outEKifRiRBWEkxlfyPJOlKCNUG2WvdXw3m9MaSlqe8qECwjpMpdF1lOwmZ5wvrclDfb9xnvvjlg82SsjilFzNi+O2J89SU8+dBl2JcHl5y+fWSdClRCkS4hnjyZAsYLH+5jR6aShHaKzz+MU5WH/Yz02EuJFBt2iefaI7WETJHaFh442uWotH/6JR9BopgrOzGqZFzkpqTqDzDzF71C/BgR/fOlOxLRpwB8CgCG115VG/KK+X7l+lyoowYqU0BTAkmaNxKARFnkAjCSgDsFHmmIgptdyAcwRn1d/mGMMNgBq0DE+PF3rLfGjTdFrPezhhAGGNLzKQLH6S1zySNBZdntAAUA/obkvjIwCD2e8v1kewNIcou+pj9HSNcJdJ3oOpHISXRi9+C10J0fyMmVXhemk24/mmxx+hwo+5bC81IEK4bsyI85FaeDmRvEM5d7X7HY3zfYvW9LMGdgeDJifLBNev09edVguPbkKEuEF0IVwn33CPbCeZ7MCNABkcAJmTJOlRIPVSBWrDxSkUyhNraUvv4p0jRHpNS2Nd6o1SSqVn1GKW5Eqpj5i/73y0T01wF8C4AvEdFHvfXxUQBfbuz7aQCfBoCL3/pmhgA5ICB5OFSts24diAAgMlRAJQEW9YQLi1xCHNKjaQI8pJv54C3dFoC48XjcNhsayBEDM64Jq0DEggDyIKSSdGshD9d/KgLHKS1zDsBCPqTBAVzc/WF9uAMNK33K8s6Aplg/UrpOdJ3oOpHKqXTiwQffZKuIQ7iP8J8nKmHntB5YESZeR7ZI2K86duGxyokqPCHZABj81DAjwjhRraEXchmuGPaBy3karjghYebAoP2I8aVdIJx248jU9n3rryEek8mRvMMF4XDfEarxwl27S2ZHQtRDqG+ryNQGaviDjFjpnnzaiJkjU0vLZLFFpE5NorLtS1Lkjk5UJ6IHRPRQlgH8+wB+BsCPAvhOX+07AfyNZQecWa+UTSXfynZdlifkStnUOpBO9Bob1rgthEMoWtoCGkvAw43tY8s/uD+puyWLnQcoCZMEgPJegDxRN1yH8hZI3kvIg9HAQGmCrr7nIixEaVJunmsT3oEsq3dWJP8lSsSTipJaK+k6F9s53b7guzpWuk7EbV0nuk4AJ9YJUh4S7zGxmefEJYz7vw2pJHIKYy3J37iNA1pKXTburxpHz5LjXd5UHBSU/NhNZBGGHgjLnoRJjtR4Ecd9qk+OrAiEdcRKponR9c2ewYMBb1wZG+BwaRz5qoQL7ZZw/cDg+mXC/iUX9nPHd6TKHRR+rCpgvIBLXL8A7AUwXjDGHWPcwf8xrP9j+dtwJJl6PCvK/zh4s5xzmMsyqUdTieaxbtoBpGy3mCn5C4SPAVhvdai/oj+EbJuQm3iqPgLgr7sZ17EB8JeZ+f8kon8M4EeI6LsA/DKAP7r4iItAZNradj1t0l3W5pTk1rkrS7/Q3KI3CZBE0GiBx5Zib6LYlTyKVdcxgmCYgmW9Y+uscsakdR6Owf7aGiGP0V/zTS1zCz+ZLJUjSrO0KkC0uhEVpBnuUPvBH0fc+6st7Wy9sORvLl0nkmN0neg6cTqdYKCcfoWTn+S6iWvlqbdK6sl6IDMM90x1HeG1vk7cr0F2iaKzxj9Xyb2yG4C3JuRcmb0nYJnHSietS57U5jGHbXRwpMJNb+PIGgBsrjjZHwAO9wj7e4TxHuFwzxGn4Row1yi8U3brSJOMpM4bTodD8B4sDJx6pcKL0s8h/9VEiNt1oLhtw6BY5YnKq855oVr5VLVjKTmaVDHzLwL4XZXy3wDwB489biGVh3aq5NzpnBGbWOPAPIgAaXKt9GrSeSI18DAA8kFwB7ipLZy4uqN88TMgIoDgdwXIwE3XkYY8Rh4ccPjjGuIQlliTRyJgoUeU1vOfjf7jbIY7kLXt0gAlFk4KJHqHIodEyubAowYmN5CuE10nuk6kclKdID+opC5qgRsnt5beF6f10jolkdLLJDFYFhIk5e7AKWGLB5fnGF4HuYm2xSN0uDAgdiE7fS7tlZIhFkAIc/eRdZ4qGQGdjSdU8h0Yl4j+6MMG5uC9T/fctW8ewSWjExLvlN0JoWJHqKRXX554PnB8sNrrE24aJZEK5ZV1ZN9bTrhwAxK1lkDVTjN96iDnM6L6lMw8SKD+sOes72qdAiTS7SWITF+bA404cnTsgi5gItcWp9+Ix3a/FgImHEeKroDI4JV98HkVI3yPJfa9pMiVjeyu61SWuSyLZS7PRYBDknNZ5Zaw2saEFBwmGv5F64hlYhmuAZPnQrpOdJ1Ysw68EDqRT8GScZdCQpJ2caC8nioqiJYqk3dUJVzkRloXwqWGVwhl+gI00SLGk1cN2BiYA7v8phFJDzyZtsZuyHm2mEEjgzd+KISNr7N3F2a33jPlQ3zv/naL7dsGu3cBc+VJm6HgdRovfJhvxy5vytsryXAIAwN6SAS7kEwp4pRvn/JGnYJEnYxALey18XyQqoosscRrZYvqoNynlpAb1ithj1zC3GYZeAhwDERlqAPAyOyAxTd6O7JV63xgCuP02NCLyUIGHBlgYNnlu9yWZS7PKk/OhQILHe4ogARQQOKVsQUsUmetJT4DQGeII4ul60TXiRdaJwjN8F9L8nGigDpZ5OZKWZaM86SJFYCQDB+Ili8TkmUZ5gCYUREtdse89xs29MqT6WGI/ThTnlyZAwdi6cga43Bvg3Hr6m6unPdrvBBCRTjcdx6oy1833kMlHjKKZOrS5UsVY01tPJGSP8ARqbVkKiNWBZFSr2mKSAUStYRAzSah107QIE4rlOA8SVXlvtaMaDonc5Z0q06RZ6KAQhJy47bpc5hQL4LHAILxX5uVSTJJrHJWjV20zi3InUu0jjcY4Rr5Ec40GnnAAHsrlnkNOCTEobdJjycgBQ63nipDvr0Eijg2T6yD6fWaLKlzLtJ1ouvEHdcJz09SWeY8KI4zJ1PqkBC1Sr1kX0U8iFkRrDjSudkDwx6BZG0f29BLjwmwA4GMI1RkEQbPdF4qg8M94/OjGHZwieqHey6Harx017N9z20fLwhXrxGGJz7MdwGMlxzm7OPBhfwgZGpjHZliAkb/N0WmwojnnJbXPFI1b1RO6CdI1I0I1FS+1NLyhpwHqTpCMYC6ZX6s1BJxF+1XyTERK11PLFsTo363ZByAwOV6DEQOJDxg7IgSEBm8v9oov7UFYyRbhDwGnM4yB1w4yHqT0RBXk3NlGwNJyKMMgaQ5JMjAZA4YapY2I/ukiIEccGpywu/pxtJ1ousEuk4U8pSuZ2Gkp77vwoMSA3vv2aIDY/ceuyR0+RZGud30iKGX4cgYLwc38bKfhHl/zw+TcM8PrbB3PfvIcvCAXb3O2L5DsDvGeOGIFG8A3noitbWgjXXfkCXwwZMpS+V3pwlUMoWMq0hqeRGJshmJmiNQOYGdeN7VOq2y1v5z++BcSNULLmKt5yGToh4MBiIMPpPEMPuxc0SLIojY7Hf0YYcdx+k6dMjjGMvcYvBgI2DhJpMdPVhIuEPuzSaWdmpJ19anLG9Qaa3X6pThj0pZJtV8kTO20F9E6TrRdWK1UPYrTLFgjJnk9ebq5/vqc649vi4z5W7ye3joPEj3vsIw19ErFZLhlYw797t57LxUgE92v6fGnaI4TIJLPqeQdM6G8fjje5h3Nm4IhK0FtgwyDDNYkAHsSOCDAQ6eTGk7SRMnA+ThvTjkga+eh/hyEpURKEB940sI1C3nSK2VTqqekiT5JnqZyFnRXgyMy8cQXzAbl7jrQWQI9QBAJ+lKt3APWOxCHaOzxVdZ5paHBCDE8hbgkHAH4EFRWemUAQuhAh7IDJ4cTGqypIHPAWbxfgvqdDm5dJ3oOrFaqLK8FBun6k+Rp9bxhTDlxz3yGtkA+wcEPADMwREZyb2SPCl74cbX4g1w9arB5Vctxq0bIkHm7SMLkB9z6vAAYaJnF+5jjA9HPPzg+3j/8hIY3bAdZBjGuG/O7g14b7x3KrsX8UiJV8pw4Y0qSJR/VrNeqGqoTy83CFZtPa+/RJZ87wvqdFJ1prKlAXtGAAYJaVgCwOx/Y8hjhMsjGcmGBN1hpWXuMkQcAIzwIKPyPgaKwAFf1yrNyS3zch3J2DxAwzqfyw2plNUs7efC+u6yWLpOqIfRdaKUmgeq5TXK6+XPoLZfLhlpWETuWteYrTMB2MhAodIzj/24US6J/MlHLC5/bYPLr/ohEAzCQKP2Ig5yai/Ykan7ziO1fXiFh5dXGEeDJ493GDZ+Hs3DAHs1AHvjrkmuS7xRBoFQaW/UIhKlCdQUeVpDnJaSprlve8G3T8W52ufupOopiUyvYUHJVBsju4EZ46SuFsbb3hpEXKKsC39MhTxkstnRdeFYZJlb383cWecmsc6BCCQhZwSpJS5lOZiMMx99EcpYKseCQGW/Ksh0eSrSdaKUrhMTooF+qs5UGU/Ua9Wdqj937hX1zMGF6/QQB4f7jkxZNVI5DwzsLIZ7Iw6vXuO9X76Hy99wie+8AQ7eI2W3wHiPwZcj6MLCDM7ttH+0w5MHG/zLH/oyfv4rH8b19YDD1Qb8ZEjzphSJgnHeLEeouEoM81BelUAd45Va+Rzn3kNJkI47Tks6qboFcY2wY/+WDfYYMNABlilJ/pUu4q4ex7nHiAG2GMh4YJHJXWWy13rII1jmIOww4smMZe7G7HHTh4wCFjkwFJZ1uj3WUesZmEideUCphDxq4Ysl8qJY3y+IdJ3oOtGlIuxe5+4d4OpVP7TBDjg8YIz3HJnS7423DLocsbk4YLs74OG9K+DV9/CVz3wI2/fI9eTbMuwFg3c2kCny+VxmGHFxb49X7z3Bk3GLh/ee4Nd+/XUX6hMxCEMo0GAdmfKbA+lnqGleFIHS5KlFqIDp73fOQ+Vllhyt+dbX6sVE/fMgVUcqujdQTyKjz5OYEtdVO9sPJuntNLKzcN00GhFI8nO5cXngQwoAwBlIGN+13HgwoWbIw7KEIaJlLnkkIzuLXFvmMlWHAxSnbTrcISI9nsK6BydkdezMx10Dj0U5I0vlSLA4W4sc6DoBoOvEDeRF1IkXRSRaJuTET1i8f4Vx9SE3JpW99JVGgA7keuhdWkeodgdstyMeXF7j6156Gx+6fA+/+Huu8Lmf/yhoJJd4vrMwGwszMMhYbDYWw2BBxDgcBvx/X/oALi6v3fdGjqCIB4w2jkiRESbn58Cz5P40gaqRKb9P676LZUwQpBOE7m5U/wg5D1KVC6NwL8o0DmvFVoChVlbbTw/em4OHC1cYuL5BHhTgAMNmoCL1gRjScADjGvzRI6EAgCu32HOsb0LvpRjyGOE9tAAKy1zlkUivJ4EyByyxK8rgAyxaJIdkrSyxvrscIV0nuk50eT4lVyvy4b2tzKvnwnWH1w64fPUJnvzGPQzvG0eQBsZ4zwIXFsPliGEzYrc74MGFI1QfuHiElzdP8G9+8Bfx2sUj/NNffhNsCcNmxDAwhsHCGGeSWCZcX29w/XgLfrzB49/YYfPhx9i9coXDfsDgSRQD4JFgr4cY0rMZicrvLRDFvJGaeRZLt1Xr355O3cSwOE9SlUvFfKpZ5EUoIVuXMpPVcePJqDrIRnLOwGRU2y1LyMCGbSF8AZe/YWAwMDtw8dckljnYJavC55E4izuO+gzf80mHPEaMboycGctc8khGFZoYyCbrN5ElFvmN5Ngwx12QrhNdJ7qcp9RIlJ9bL05Q7BPO/a/dAua9Afv7G2zeGUDW50Xdc96pYecI1XY74qXLK3zo/vt4bfcYDzZXuD9cY0sjvunhWzh8zODnv/JhAAhkigEcRoPrqy0OVwNwNQDE4HuMD7/6Hj58/138zK9+HfZPNsDBxFCgVSSp8ETV7nv+u3zqXtBn4HV9PkhVJjULPS+rgoe3WHWdBExAsCo8McLAMAfwGNmEhFS3vwmDGjqwkFwMZ68PNLpcEm8lm4plDhVekJAH2IbG2c2J5jwAgOtuvsQy73K3pOtE14kuz0gqJCoQqcEPaTB4IhUmJXZl8N6h4colmttfuwA+9hj7d3aAYZj7BwyDxbCxuNztcf/iGq9dPsYHLt7Hg80VXtk8xtbrlAXhX3n4JVzbDX7lN18FAxitwX4/4LAfwAcDGhjm4d4NMmss3n1ygUfXW+y/dgHam5L09FDwajl/UrUgyJ9b6M5C1esEqoBFankr4GAfxsiAxKhEWwtOgEZCG3IcscjRsMwdAJAKbcSQhwMR+HBPtNAFRFyAJVrmowc52xXgbkjXia4TXZ6dtEiUJ1IwcWyoMH+eLDvnqptqhgk0umlrABmkk3GxHfHhj38Fv/HOA4yHAZvtiHsX13iw2+PeZo9Xd4/xYLjGS8NVIFSAM3D2POCN+7+J9/c7fOnthzjsB1jrttPGTSRuDwTeD8C1wYHv4fIj7yfX0eVmcj6kSuK1lJXl6xmg5NZ4sZ6fhimxxi0TLMXu3G5brJ/UBSV5JDqHxFniFiD/O2OZu5CGi1NLyENAZEduoEM91YW20AE3J5rtZsSLLV0nuk50efYyQaKS0J5ByJdyHir2BEqRKf+p00hBZdnA9di7tMDliN32gNcuH+OViyf4tXcfYrcZcbk5YDAWD3dP8Or2Me55UiU9Svc84MpuYJlwsAM+ev8dvH+9w6+/84rrUuunmSE/f5/xBIq3XV9OLedDqnKZBRPy1vg0WCALfyTAoSxvy5T0dnJhD/LWrgcOHfIgl8bqLHQHHAN8ryTULXOwf9w0hm7eAHDNBjuy4R6v2XcPzyz0EePt5mt0OW/pOtF1osvti1YaT5xkOSdIgUiJJ8pPShzKKT2GzNsH+P03fgyqSwtz7wBjGI8f7/AbF/fxOz/wRXzDy7+Gz7z9UTw5bHExHPCBi/exMSNeGq4AIIT9hFDtecDjcQsAeP3eI7z/yg6P3r4H7Am0J1A+f1+Xk8t5kqrC8k5DGQWWZGAi61BgYdQxCyApgCPLIQm9mSR/ZAzhjhCCYBuSc90+0TKXOcSk7sCMawJ2bMONaBAZgRjC8MuGOVjiRnpFdbk70nWi60SX25GaJ8r/am9UEuJTHioJ5yXrchxS52B1XOO8RPbCeaeGi9FNd7SxuHdxjZd3V3h58wRbGvGvvfpFfObtj+LV3WNcmAMuzAEPhyd4d7wEADyxW1gmbIzFb17d92F8xsgGD+9d4cHlNd5+7x72VxtYmc/vQDDX2WSEXU4iZ0WqiAkss84D/gvmBC0KMFHrzBTGmwFi3siQhTUkMVfnkBTAoXJILBNGPzaOJWehw1vmCLkk1ueRRMtcRmrW4YwBAh6ogogbq4dTlPTJt5Zc8u3Icv1xeex48kJK14muE11OLFOeqCy0l3inpFyF+9jEEB/ECST7MpCMIhI8WgzeMvjCgi5GmK0bQ2q3O+Clyytcbg64Gjf4f7/2Bn7na7+KLY34/a9/Hh/bfRU/++jr8MHtu3hinTfqkd3hym5w31zjym7w7uECOzP6MKDB4+stHl5e4bWHj/DOcOkS1r1n11oXEhwGm7QnXW4m50OqVoY22INLBI92Yi4j7flUyyGRZQccLodEwh2g2OMpt8wHGss8Eg8mzsLmMPAhYHDNwA5ogojrjs4ecBx4umAJApBoGeG2W7h8lnzU5i7PsXSdCNfSdaLLjSQnUv6X8+XMKwWknqnooWqQKZHKRMQ8sBtOwQ/OSTvXq08S0e/v9tgO7msm4kCsfu/rv4yP7b6KD23ewe9/+Bg7GvFT7399yKMa4HrG/ub+PiwTLswB7x4usB8HHEaXGXm5OWC8vMY74yXswbi5yAn4xBtfxmAs/vnb92CfDHEYBfGs6edWCbPfeIiEF9DwOR9SNSU69MHpehHm8GUJ0EwAh1WWt02AwwC+h9KcZZ7nkbg8kYMDHFBYlzF2WiAiCbjhRgghpGHgrPKa7NmdRwYYtOG3dO+OlbJj5dZzWXquTFu6TnSd6NIWFWoLv0KSat4o/YtIllIypUJ8sk2fQ58Xen83HhVv/Px9OwuzHTFsLHYXB9zb7XFvu8dgNBNz+nttB3z+0QfwO+//Cgay+Kbtr+E37QX+rZd+Hn/33W+AZcL94YDH4xZfevIQL2+fwBDjYA2uR9d9ZD8OuNwc8GB3jevDgMfWeafZEr7y6D52mxEXD64xXhKsNW7gz4NxyqUH/4Rub+R+tYu88Q7myoo6p/vGn8XsAOdFquRl0cS6srRpCkhUDollcjaxyhGR7uRiecdlAY7UMpfk3JplnueRDBTH6HEXKgm50yAywI+8LH1bPTAMfj6DVvdwAQ8Lwp6Nn8PM/Vk2s6Ax4nSg0uXE0nWi60SXeamF9DISlSxXiFQS3lOEiv3cdwWZAlI91CKJ7ANHMrVhYGthts47dXGxx/3dHrthxGDiOGzhEMS4t9nj3rDH//3ot+E/euWf4pJGPOEtRhC+5cEv4p+8/3FszYhffPRBHKzBvWHvQn884Grv4P1qv8Frl49hifDw8gpX1xvAEt788G9iQxZf+M1XcPX+DmbrR1/fWfB2BFvj2gyLMEVNPj0NM9cJVv5M1HJ1Shpuriz3ZlXqreZnJyB050OqcqBAaomneSLp8lwOSZ5PYogTC92yCdZvyzIfiQG2s3kkIDedRWIpLQARN30GJ8eK3gOO4Y+K1MBjzwP2vPHhDxN+3fFMMW1GPjhiLWSyJIxSs9RrZSeb4+wG8iysmFXSdaLrxFOWs9cJkZonyv/WQnr5+lTYL3iXMoI165lCPA77QT7hSRUGBnlCtd0dcLnb43K3x4UnVCJhqBBiXAwHvO4H+Rxg8S+uP4z3N2/ji4fX8KHNO9jSAd/68LP4lf0H8A8e/w482F7jwnutrg4b5Z11urAzI5gJb7z+NgDgYjjgK48e4OpqC7MbwSPhMA4gwzCGQcQwgw1z/Vl2D4ut/1bzeQBFx4O7L3tGvpyXkKZK2eQEyrqtnJK573tKARbqxvmQKi/EBG6ENWo5JKkVzimQ+JwPVsdMgSNdnrPMa3kkMocYZLwdNtgDETj0cgNEBmnkyWIn9SBzp3H4WFqhhdFfXw089MBwYwCYWKaPWYCIKFJSZoo6uRybw3I0qDwvQHCkdJ1A14nVOx6323MjDRIFpCRIr1fJVkGqJsiUPm+DTMH4EdMHdhMUD85DRRuLYWuxu9jjYutCfltji5AfgECoXvZjUr00XOHDu3cAAD/x/ifwtcN9fPODX8Zv2byNPQbsecC//MqX8U0PfhW/8PgjGJnwzvWFS0SHS0h/9+oCv/Xlr+HB5hqvb9/H//fodfzyO6/ha28/iI9y4GCQ2ZEAIhB7Q4aiQQNDvm0RUuW/08R7VSNWjfLwTNW3XvnsudW7d4JQFUSsVm+prlBjOZNZUkVEfxHAfwjgy8z8r/qy1wH87wA+DuDzAP4TZv6a3/a9AL4LLl/0v2Dmv73ogpmQ92qK25ABiVqvgAonoIIi3BEGN+TY+wkLLPNaHgkESIAQ8kCSSyIPEk0QcSENByTILPQRhIEZc4MatsDDwuCah8Iid5Z6LLMZ0KTHzkCjAjS5lMBTs8wnb8lXqmnXkv0Wlh0hXSe6TnSdSOWp6UROnLLlKnnSy9Raro8ttYRMBc+WJ1OQXxMJ1WY74uJyj93mgMvtoUmoAGBrRry0vcLru0e4N+zx2vZ9vD68hye8wyO7g2XCZx69gSeXW7y+eQ9f2r+Cj+zewRvbr+GN7dfwq/vX8Ll3PoTLlw/Jcf/1Vz+PK7vFE7vFwQ74zXfuu1tQ3hnyHusQ9hsJIxHIsLtVT7CI3DqkLqTdYdcEsXp40ubo51f1YnFali8T6t/+BFFiqnzgtXeoV481ZLwsSRz4QQDflpV9D4AfZ+ZPAPhxvw4i+kYAnwTwTX6f/4mI1s1U6plvsi6L2bakmrxIRMsbcD2eZJnhltPtFNyjoQwUG9qkzFm08isDr1nfWO95cPuBsOcN9jyEerKst4/sGve9678R6u15g2sMeMLbUP8aZvJvDjxaFvno93XPwgRwcPXj55GHR3Rd/Xy0lOuV151/wEcCRk13loQyjgx3/CC6TnSd6Dqh5Qdx2zrhvUFhCphknKj5v6SeJz7We5RkFHQ3Hx9iREuTqRzgfZhPJkSW3KngodpamN2I3eUel/eufUL6IeRQ1WRrRtzfXOOV7RM8GK7wyvAYrw/vY0cj3hsvcWU3uDAH3B+u8QuPP4K//dV/DX/vy5/AK5tHuE9XuOYBb12/it/60tfw8YdfxccffhVv+uWv7R/gtc37eOvJK/jqk/u4vHcNqlwHCWkaGOQ/dx4J1rpEdmsN2JpgxJFxuVjuz+1nBobZ+GMM1v/6ZV8Hxj8nA/+ry7Lt4Y/TP8r+9Ptp/enjVba7Dgn+j+p/Uwo466li5r9PRB/Pir8dwLf65R8C8PcA/Blf/sPMfAXgl4jocwC+BcBPzJ1HekunFkG0vLXFrruO58sS7mBfTzNv34u0apkvzSMBbGg4rX9DEvKAt2a3wGLrfAihkgMKC93XHxYy5znwcKBVt8h1WCQB0yzcMQcSOSDkn14tnFGCSOXmigPV6ix4Tiew1rtOdJ3oOpFVf0o6EfKegNJbhcq23PNEQAjxkTqe7Ad1nBp2JsCLOI/f4H89EaANh/yp3faAC0+mtsOYzLepZUMWl8MBL22v8PLGhf1e2TzCpdnjCW/xlf1LGGBxf7jGyAZXdovfuHqAi+GAt65fBQD8ypPX8dl3fgveubosjv/W7mU8emWHX796CdthxMPLK1hrcH09FN8JicebGASK0TsGQBS8eMRqCqzMgwX49geqLXOFrh0KD5xVfamjrke7bxOrkVIvk5yoVlfvk9xopU5t+wp9ODan6iPM/BYAMPNbRPRhX/4GgH+o6n3Bl01LZoUXOSThxig+OEZgyjrcoXs8SRKuBYcBD8UyN9DAIfkhqiwBDgcUlhjgMkF3LzOfVUFkAzf/mbpPvyzj6LiBFzcYYTH4gRgt+xGogXj8BbIEPK5507TIdahD6gAaSJF4K8K+GchoUOB8HTXQWAIipwGaW0rG7TrRdSKcs+sEgFPrhCZByG69Rp6QrTfypYAVZIrg8qbE+zVwGLMqEqqYP7XbjLjYHrAxdpZQXWwOIez38uYJXtk8wgPjpqP59cND7HnAvWEPA8Z74xbvjzsc2OCN+29jSyM+9+gj+MX3PoC33nkZ14dSP96/2sEQ44OX74Wy8fIK7/IF9vs6sQoqSghtCgA/zlUkV0BJsMIxfJskBAvsaJoOCXK6W/L8E6KlSU5OoHSbCLWujzdHuvT+WpbZcABOn6heO3X10onoUwA+BQDDa6+lgCB76Qfo15l9FRLW7N+3fzFEshxf5Jxl7hrsumUegQO+a7e3RH2CLhaBSGmdh27hrPYhH1KAxaByS5aK5IgsBY8pi1xCOkAa5tBAcmOLfdaCr9xkVlYNcRT7HGmtn0a6TnSdaK53nZg+m9aJzcuvyWcUjqK9VXNeqXx7ctK1ZCoM/smujg/56fGndpsDdpsxGTKhJhtyhOv+5jokpt8frnDfXGNLI963F3j7cB8X5oAtjbiyGzy2Ozwet7gc9vjghSNJeza4GjfJN6Y7rlgmXI0O8l/aXOFgDeyWsB/d0AmHg6l+F4nXSsiVf3jiuQqecHm0CwiWJkGBZElTp8hVGMBYb2sRKE229OPO6+f7JDfM09/+jO4cS6q+REQf9dbHRwF82Zd/AcCbqt7HAHyxel3MnwbwaQC4ePPNeAv6zQhFzsMdyvKWF6wtcnf8FZY5HLBoyxyAT7rVwJGGPPYYsMW4EERS63yEAxVnletlVtuclb5GxCpfAh4up4UmLXLJrQllCZBML8s+0xZ7ChQ6Dyh+D5kFdSLrfFGd5dJ1outE14lUTqoT9z76Jre9ULFsFZkKJ8rWczLlbAk1kjr7nBwO4T5JSN/tIqGShPQWoTLE2A4uMf2V7WO8un2Eh8MT3DfXuKRrjEz4yv4hLBO2ZoQF4b3xAlfjBgdr8NruMS7NHns74NpusB+HcCuaUMkt7kdX796wx0vbK1gmPNj5PCneYhwxSawAhJBgMOICufL7yvNnCfFNECz96GV/d/HFeym8WbXtmjzVSFOLcOX33CJccoyJJmhJonpNfhTAd/rl7wTwN1T5J4nogoi+HsAnAPzkoiOyelicNRKqPOmSqZ+PBw8NIOJu1Mm6kpCb5EcgLRfXfUzWjQm6ruE1oY5LxI2N7Z6HkKybJ+pKo62TdV1j7/aPibwxB8SVrftbCx4CHFMWuU7CTZ9frJsnO+vv0WbAkFvcq0MflfWi3VpgwRd1KlUWSteJrhNhe9cJACfWCSaEPCaddJ4mqnOW76T2oURF6iJkSo41cJLYHnv2Ifbu8wnp290BFxf7dYTKjLgc9ni4ucKr28e4b649odpjIMbb4wM8sjtszQgDxqPRzfV3ZTfYmREvbx7DgN2UNaP7tuW7oey8zIS9NXhycPMGXg573N9c4/72GvcvrrHdHmBM7aNBcjwCIInhRM5+0p4+ZoSBQuV5h3C3/pNjGOvGxJKcLGKXzG5corwcn0y6DZ7wxuP46/F1Q/J6SGqP77dIUq8lvVNlP03gG7JkSIW/AuBbAXyQiL4A4L8G8N8C+BEi+i4Avwzgj7qXxp8hoh8B8LMADgC+m5nzqbnaIgAR8kMQWWSwyKFCGq5cEuluYpk7419b4+2QB+DG3dn6reusc+Ptepfc6yactbCwMLAYMSRW+jHiAMssBo89D9Eqh0lAJLfICyCpLGsgqeWOhOVsW+xy01gvDjCzHfOgUq0zI10nuk4AXSeS+k9JJ3ROVeGV8mWLPFPFDejjsyJWfj33Tvl8Ktq4+fu2uwM2mxEXmxG7zcFNAr6AUN3fXOOhT0x/ODzB/eEKF2aPLR3wxG7x1cMDV98Tp8fjDgc74MAGDzdXbgR1OMPlybjBaKf9JMyEa+sI2MVwwOVwgGUD3hJGP4L6fj+4capqj0l7oICQ55l7rsIz5xgaRLZvcgblwXKrFNMVgBB2DO+JxXum3p20k/7kMWxI4RyJ4UJqm76g4pWtU4Ylvf++o7HpDzbqfx+A71t1FUCIZpBY1iQPNgKHTsL1byw8VIavD16VRwJE4JgLeeztgK0ZA4hIku5SEHHHjyNEu/tyeSUjOyAZyCZlx0gMR6wDD/EwaO9EyyKX7VImQFB4PCoWew04ks92gQWuZZE1zun2fJ81INJ1outE14ns8E9JJyS2kpApTZ6OJFOOREUypdeT3CnVw89sLYbBhfy2m9HlT60gVOIperh5gpc3T0IO1Y4OsDD46vgAex58LiHweNy6vCm7wYasG2mdLEY2nmgN4XvSkudVjdbgwAYXADZmxG44wILwgK8DsToAbvDQyvGSgYblGfs2KBAo9fBlXehU7HATj10QLO+xgiZw6leTrBAW1LlQYmCGa5M66XkLolVeTWxrF8r5jKg+YYUnTzWxuHG0ZW7gXq5lgvGAoruTA5LkanGwBhtjCxCJy8tBBGQx8uCABBQs9AEpkFg4RZLeTktFEnNLC3saPDQwhIRdxLBPyyLPASOElQAFKnE5fzd57kiyPGWB5yCQbS9AA+X6Eov9mUrXia4TXSeiBA9UI18K8XfRZev982lpAsHimCSjeveZDWPYjNh4MrXduOEStiqUVRNDjIEsdmbE/c0eL2+u8Nr2Ee4PV9jSiEvaw8DiXXsPj8YLDGRhwCE5/WAdcXqwcWFCwCeg2w2uxyELJ1MIp8k6M+HABtejCx9uzYid8ZOnbwjjzvhvcwvmoUI41OOb8lrpdzBDrghRF3TuliuIF6BJVtWL5Q+2hmQh399vX0S0GnI+pAqpRZ4CCuLCDS1zwJVZXzdMx4F2yAOwYdsxIGJBMHBd0g0bZ3kEy1wApQ4kx0zsqvM/1lrjI5zV47alOTRTFnlrHB9tgWtQSco5WkR62RXkwFACiZbiW88srbMGjYp0neg60XUiCqscl5OSKfmsJbRH+lychPvMxmKzGbHZWOw2B2wGu4pQXQ4H3N9c4+XtE7y8eRw8VFs6YCCLJ7zD24f7QUf2POCR3WFvB+zZYGNG3BuunY7Bfc8HyQFUZL0mlhE8VQc22GKE8STPMuH+lhJD4HAYXH5U61GqNiQ2S8eRK7+SHLsgWIqIhXMqguX2U8dZQrLk2FqnChJVI1p1OQ9SxVDAEPHBLfsh8EnX9YsrLHPAgYYDCmcRz4U8ojVucLBwy0eACGACYIxAABLr50Zzc5mVQGJXjMWjRQ9iGEBjsTVOCjxMBUjKZXeumNgs357bHpfD61bvJglzaEs6BwogXU9AJrfQ87ozwFGz4J+1dJ3oOpG92zuvE0A6xhTi72oypfKmNMFKyRTScN/GwgzOQ7XdOs/UbuOGS1jsoRpGXG72eODDfi95D5X8AXAjp3OE5iu7wbXdBOJ+b9jjwhwC4TrYcigFQHmS1DXpzicHazAal09pibEbRlg2uL+9xt6HAdnrwBpiFd6HCgkm76hBroBIsAJZUudY7MHy26vH19edfTWM7PiFcdJ+BlrOg1ShZpGLFa7KjrDMAZ+DAgqhjcUhD0RL/RgQGfxcaAZxzJ9BWeMQ66IBJDcRAQy9XAMP3VMrtcbF+olAoi1yWW5Z5LnVzfKrLKlWt/EkzBG+BbWulqeAg2a2F0pzZiDSdaLrRNcJJYoEiRxNpgiKoHHbO2Wcd4oGxjD4pPTtwXmmBouNcX+LCdWwx4PhOiSmX5q9T0wfMZDF+/YCj+zO7QcOPf3ES2XIhvGqRPbsEtclH2pOJK/KEasBw3BwwzWMhN3g5gscd9eBgAFYRKyA+L3mXqsl5MrtV3qvNFkidbxYIX7UyRXWCFbDi+VWJ0hWcaVtORtSBcA9ePZxUiGmfKRlDnlgwqJlwLH5kIdY44aVZb4QRAzZ0AMqWuWAC18gAIm4docGkEjWiKHl+SMCGACSUaF1d/YWeOTWeAIenIFHxTpvWeStZFwJadRDG6SW4zsPQwoUwJBa8HPhDMrqnhVw5NJ1outE14koN/VMSZnq5Vd4p/Jw38Awgwv3bSR36lhCtbnGw+0TvLS5imQKFjs64Jo3eG+8xJ6HEBoXUiXf2YU5uP38mFWWCQfvrdJhu/AIsusSMm/h9rMYMbLzVm3MCGvJhRc3+2gM+H3niJWcryBWgGq8EFIQUrcWEpITvEKJxwnJ/TUJlt6u9UZOlxA2ve8KkjXxGM6CVBEAwYebWuaxSnxza0MeNwER6Vo+gGCYg1U+sO8NosIfIexRARLAgcfI68IdCYioJNpynKBWvkgJHgcNIgo8tCtZ/01Z5PIe9NhuwQqfsJirVnq2D7Xq+fUcOIrznRGQdJ3oOtF1opSbkqnopVKESkgUEJPR/eS/jlD5HKrBYiMhP+JFhMrA1duZAx5srvFguHYjpvs8KkMWWxoxssEju8MTuw37P7FbXFk3gXjwUg1uyhsRGZ+q1fOvJvq7PFiDDRkM5O+JLGAAy45Yae+XtAlriRUgbVAkQcFrtZRcISdDCwmWeLA0qdP7KtJWnC+/r0lljHIWpAqAAokjLXMPFJSDCy0PeeiBUnSS7jEg4qxyVxc+F0SAJIAK6kBivaVybPdxDRqyLpZ4sMCPAI+DAh4d4tAWeQSVuAwggMpcMu5smENtEz2qggHKeq3tzW3PWrpOdJ3oOhFk9jKWkCnxThkP0jXv1ODCfca4/KlhsNgMFrvNiM2wjlBtvYfqcjjg3rDHg82VSkwf07DfeJEkp1/ZjR9U173vCzNip0J/0btqEo/SEkmIFRsYthjIDW1imDxxO+ByY7wXzH2Ph8OwmljJqwnXdwy5EqPQ76+q+P0pnDdIg2C589brNb1YWK4GZ0OqxCJfbZn7Xwa7Z0Z6f4omONohj9HXyUEE1gAeJJaCiCGLAexyUfxHOnhF0eEPy9QEEgC+/vGSWONynzcAjyLkUQGP3CqftMh1I5AAA8X3noCFUtBa/VAv1lmcO1IDmTOQrhNdJ9xy14lJEWKEOqHKQ31V75TyUJmBYYzLnxp8zz7xUq0lVDvjBth8sLnCg+Eq5FFJ2G+AxTVvXO++LOx3sAMOPpcKALYmkjAAkHC1qxe/RS15orr+FuM3bGE5equsv35DjMvhAN7qb9Yx0qXESq5BXhOQkatAhuR6KxUXkCs5ftV7BbWPnDfbr1aPi3rT9ytyNqQqoZ68wDJnhRQcy7Rl7uYokhfSDnlImaznIGK9Ii0BkdQid+uWEYFEwh/EVSCxPMAQp9dwpAhQyHJqlZtV4KHX8x5SOXBoizy3xsv8Eagy9R0E8tAo0+VAO8yRAYcGlhxkzg5Euk50neg60ZYZMrXIO0VIwn0khGpgF/IzNhCqrR/Qcy2hujfs/STJ17gwe+cNojg46CO7c2E+8VLZ6KUS2ZLFhsZwfJl54OA7WbQS1POcKiAl+gc22IRvmCreKuBiE799ALjGBoDLh0TjvPk1LPFaubp+dSW5yqr5YywkWFlYsVoPJclqyfmQKiyzzElCFrJN1RHLnP3TTxJ0UYY8QuIrnGU+mDqILAl7GK9oprDI2VvYEVgsD00gceeLL28IV7lM9Bg+ogQlkJwGPKZCHJDXokBi0iJnih8xVyzyDFymwhzUqoe0XrJ8puDRdaLrhCvrOhGE4u8iMqW9U5pQVcJ9ZnDeqcEnom/VkAlzhApw3+ngc6ich+oaD/x0MkXYDxZP2OVNjWySsN/Be6DEu7oxcb+kY4V4qSrESrxU8Tdu057VuJx6q8BOhzcUE9dTI8DAWhxNrIBpchWudym58sfRb+dUBIt5qZ/qXEhVrsA1yzyxvmOdxDKHNDws/yCPJlrnDGsNjHFpuEtBpGadG78uOSLp9TugSIBEAUsNSEaGt8w9kBDjmHF5tAs4nwg2GW/HK9Ox4HGwPoBTscDHsJyCCEQpuWGRVwAh/DItCnMAGWBk4FNdRrb8rKXrRNeJrhOlVEjTKjIFVMN9RDEhXQiVDOo5EC8iVBuynlCN2A0j7g173BuuwyTJjkhxCPuNMLiy25CcLnP47XnAXoWht2Sdp8rYEPoD4lAKediPs0tMxqlC/L7SMLYN5xNvFWBCoj1wwGEwylAguP4ixxMreZ3J5bZCglIZagdK73WJ9wrq/gsv3gqCNSXnQaoQn2VumZcWesUyR/zNE3QZsi8lB5JcEiAm6a61zi0ikBi4MIa4deXPXYC2yKkJJIbi4IYCKMeKBg0AQWGk7JTgIT1EaiEO904iSMSckvg+qxZ5Bg75e04s7qx8bZgjD3mci3Sd6DrhKnWdAOBAVOdHYSWhEu+UcWQKJuZP5R6qmxCqy83eE6q9HwLhgK1PMNdhPyFUOux3sKboybcxIzZmTDy0lt0YU/kQH/oKa6E/EflOAUSPFcVzDsQAXBgQ5HRvN4zhO/dX5n/XESsgJSY5V5oMCdZ2yLZPkatkd30NKwkWTdzq2ZCqqruZSUzy1DKHruuBgOMvg0GW/Pysvow4JtmCYC1gfLKuWOlLQCS3zpNABAHaQrdM1fBHDUhGEAb/ut3AiDd7nBbxrYviyTXVLHOdB6J7NM2BRzVnBAiWuLUmscjDe1VlUxa5noJDLHLSAJOBy+IwxxQ4nYt0neg60XUilZt6pwxcuM8TK8mfGgYLYzyRMnwDD9UhEKp7wzXuD9dhXj8ZPkHCfjo5XTxUV3Zb9VJtfT4VgGTAWgkRznlPtFiuJaxHwwGwnlQhGEGGGDtzADbAweuNGAP+AS8mVsAKrxXiMWfJlS+bCg3q6i3v0xqCVZOzIlVixMpM18SoW+hAO0HXZcCimktCHPdDTNJdCyIyZo+AA7PvMSF/ykIP0rTIXZkbaBEBTE4hGjSAtmV+E/AYVXfb0W+rhTgKi9y/c/ebWeIKNFLLPZblFrkOgaThkMYyymXK/ebPWrpOdJ3oOpEIy4jqiadqglDJsvdM0cAg4/KnSBGqwfhkdL+8llBtjMVuOOBy2GNnDrhnrp2Xig5J2A9ACPvt/XhrOuynR+KXATnFu6VDf5YpkC+pv4ZYheMgkinLMssBwvlF2yWncUMjLoZDMBTE0AkhxRMQK6BCrhLSI+es7MRxPSFXUNc14b3S91L19C24t7MhVQIWBV0N1na5noQ8LDlLJFBVaXDYdf00DkQEEBgMa904GcdY50nGCKFqoesQRgh9NC1yd4Mme2lTCl0Tm+0fu5GXlnluoeRla8AjggPB2rx3E8B+36ZFDkSg8MuTFnmoD1U/Hj8pQ2pgnH2Yw0vXia4TXSeUrPVOSbjPEykJ90n+lPNOcZVQDb7n21JC5cjUmIT9Ls0el2bvcgeFGMHikb3AE5+cDjjP08G64TxaXqo89LfPiH7+jc9J4aFCNDIkYV2L8caa81iNsIP3cm1ivcPBk8QbEitg3mvl9lXESnZCtmPuudLHUO8252dASlKnQqm5nA2pAhAUnhEfmFjW1QRdvw+A9C3IJksgDw6CQTmIHGOdSyNqiF3uSbDGSwtdwh9icdeApGaRi/t1rbKIaOtbH6c2J9kpwCMm4WbA4Rt1AY8Y9ojb6la3lNe2+V959QpIUqBRDyQBqZnlc5KuE10nwnfQdSLPqVod7vOEyhgO+VMu7Jd6qIBpQiVeLBktXQiV5FDJdDIyfIIkpz/hLZ7wJnipQmI6D8lAn7mXSl+L62iRfp/N58W0iBDkBEsS1rW3Cux7A5oRO8TQIXuCxQNBBmRZS6zkWpNyuYekkAti5fZt7KjawmPJVX5tc8/zPEiVbyQInIY2BERQ/wWgQhdAQBGOaBKSdC1ARiCjtM6NiRczByIWvheSf9A1C31kSsMfxFUgkSRFIFrkAig3Ea1owQrhHDxK4GBOLfa14BEtcsBasbyjBc5WA0YGHnLNSaNOGShkFnkCAg2gycuqyxkhedbSdcIdp+uEvIL4e1d1AljunfKkSsJ9LuQXyZQxHBLSB7POQ6Xn8xNCpcmUeKh02M+Ns2ZC3hQQybwbuDMSJfFSDUpPJPQnXi0LH/7jSG5q+lEjAHnvv7juQ4BU91bp57Ih563SeYqqlvfSpgRoTm7qtXL3lu0ILCdX6ti1w4iILrfkPEgVJqxtEQEFztYtQCbLJfE5JAImDD4JiBABTH7MD2+pSMNaWOjEIfwh21tAIttPNXK0SAIiCXiYAlRq4CEzmR8DHhE0VM5I+CMFJEg/es4s8lCW/fpNulx0rWaR1yz36vIZSdeJrhP+0XedEJHvcc47lYX7JH/KJaY7wiKEauPHoDqWUMWQ377Io8p7+11lYb8pL9VGjWclYtmNryaESn+3S0QTqnBM+b4p1Y2WtyqGAV1ulR0o/GpSNI7m9ogVcDJypfdb471qydmQqtAoeIwg+RXL3ANMEfLwO6RdyuWJRTA5FkSE6UuxhUu6nbPQA2hARsAtgSTpNp5Y5Ms/wjnRoOGufx44dLx9KXgUybdWgYcHi3rPptwyT7clYFCt09gGRODQy9n26vK5SNcJt3/Xia4TQOmdCmWceqfCH0L+VPRQ+fGeDIdpZ9YSqo1xo5vvjEvadpMcj8FLpYdPEEJkmZKwn05Odx0h2l6qMvQXw3Ty7YZXtkJPWNXPw99O56K3ShMr2Q7A5VcJoVLHc9+9az9OSaz0vcYN5fHXkCtdbwm5mpOzIlXud2HIw8JNGaYop+sy7ivoJN2VIKITdeG3O9CK1veche6s7HkgcaGOaI1rF27SU2qFCFgASKyYPNQRFGoGPPIxd+ascQ0eaRKuf1UFaMiXnVrkhSUuvyjLq3W91Kzw1Irn87PIga4T6DrRdSKVqneKUA/3kfuVcN+pCNXOHLAbHKFy3qo0j0oTKjEInvAWe7txeXlQI6GzDOURvVSiG5KHpUN/rRA10CZUc3lVlgkDSWeNGAJ0E5hTyGUMzyJ8cEjyq5gJdhjDMaX2scSqdU9LvVbuOBmxkgMA68iVPseMnA2pClY4p6BQPD12jzQBE11HP6RQvg5EAErG7AlPnOKAiEDbQhclmQISnV+iLXIRsdqPlTzMoX+XAgcn244DD7YKPHw+SQIeMilnBh6lBU6KZEDtM7ctLSuWc3A5I+k60XWi64QWRaiEXGXhPk2oTCBTHPOmTkSoNJna0hhyqCTsp4+5500yhMLIBiNMnLsv81JtGwnqgBCflIhJefWJVco5+7YHX2aZglrLeYC6t0o3ARsacTAGW4xKr5IneBSxAo70WgHJeapeq/wgc+RKH3NGKc6CVGlrS0ABCaAoy9wCHHJD/P5ijWfLKYocByLkLW+iuOyua9pCnwUS77Jk8ShkQHKKcIc+Rg4asr1IWkQbPNIxd24AHjp3xL/rHCxSK5uSbyQBB0xY5FnZpEWur+cMpOtE14muExURQjUR7iOjyNQQ86U2qoffxtgYxj6CUF2YEVuy0UNFh2T4BAn76d59QAz7yRAKLS+VDPaZj00V86nUUCATetEiJVoSMgUV/kPbW5U/o53xHqrBGylDfs7TEyug4bUCsDgkKAcB1pOrhpwFqQJQsaomQh4M6FwSqS87nxJEEkq+wEIfeR5IEoscFMAEQPiITyU6X0T/rgEOVnWOAY+0V5P+i2WkQKUKCBXwqFnkNRAJvy3AQnaOc5GuEwC6TnSdUKKS0UHt/Cnp6WfIe6bEW3Ukodqo34vQ228fvVTmAD18ghgBex7CmFS6x94ICeHFRHMDTrxUOkE9jKLOVPzJ91eTJYRqgNMJkmP6EGDsMBLrN71VZvT347xVG6VTTkfcHd4WsQIq5KritXLH8sU3JVcNMXMViOgvEtGXiehnVNl/Q0S/SkQ/7f/+sNr2vUT0OSL6OSL6Q/OX4PerNhq6gUHZiFhlrdm8vm6wsmVfN4wXE5JG3bplcgPzWddQpgP3wZfFePhoDaw1GG1qwcqfNL4jU4yjq+0H/2fZzTV20j91bMkDkWs6WOOuPwGNafBIn8cy8Eh6NSXdx+Vdtt5b/u7yv9x6Vx9U9i1FK3zGIm+3sfFb7TrRdaLrRPqtPg2dICThPhrin8zfpwnVxg/qeQpCtRvi0AlbPxffhQr/tcJ+Mko6EMN+7jsbihBeCIU3rk/yqQAEQjb7yCbuUUhKi6xEo8NUjRp9vYYYGxq9R2/0uWvOO7gxNnQSMMZ5sY+JK5M3yprbmxu4ej7naW8cqEawAIRx5iZkiafqBwH8jwD+t6z8LzDzn0svkr4RwCcBfBOArwPwd4joX2LmcsCLXFhZ3+oX2sIWxfeLBPjBDL0NrcMgWG+di9LK9B3OMwDoZF0Jd8QrAOAt7TT8UbfSofYUS12YNKsXv2YE1+JRZgqQx95ZLwfFmbfEmSvLx4KHfxCtnBEhEFMWedN6V2CQAEQOFApU9O8C+UF0neg60XVCyw/iaejEkHqnQDEhXffwG8gPn3AkoRKykBMqSU4XcqUJVTXsp5LTxbMUp6NJvVS5p0tL4t0STytSD1X+jTNTk4hYhsulUsc3co2UhgDBqHqrclInz8sawhZTr/J4jxVwpNcKQC0k6I7XIEoTnqspmSVVzPz3iejj84cCAHw7gB9m5isAv0REnwPwLQB+YtHeuYIjNhwMhFwSLQno4GYgAvjpO4jdkzYpiAAUQi61vBIJf2ACSCTUocFErixJzj3iY8slyR/JymrAIeVLwEMPYngUeFgFHt5LkgCAnQaPprWe1QVaoKNISb7/jHSd6DrRdSKVp6IThMn8KfFSaUK1Hdx4/GsJ1RA8LiM2Pnlcjx0V/ibCfno+Pz0mlXh+ci9VMthn1utPPFzaiwqk+YDF41rIiCUEmJa5EKBcWz7RsiuvhwHlmNZ4kmas1xEGBglp3h6xAtJWLd0g33i6bzMkKAcDFuvDTXKq/gQR/XEAPwXgTzPz1wC8AeAfqjpf8GWzIm2mgIQDAr+RfS5J6DKeQYEHicUgwki7lpPsL1Ah56NQneTiSJ9dGuMcVOpAAqaQvCvbNJgAsafUKaSWlCvlEUDitsWWuABEAzzSHk1YDh62tLDrf2WIQ1vW1WOgPE5e9waOEJGuE10nuk6kcjqdIMQhE4BJQiXTztyUUMlYVLss1LelMeRUDUiPKWE/3dtP50Pt1ejpgPJSIY6+3gr9xb9lvf9qInmPmkyJh0ryqmI+VTrDgciQ6aQeZsH4wVGrMhrf7liwnx37NokVMEGuKvsuJlcTMh+Urcv/DOB3APhmAG8B+PPZabXU74noU0T0U0T0U+Oj98vGQP/WrDZ5h6EOVffJt0fzVDVqvpFjq3NJ5NfVYaCaVzKVWyL7SH6Jy91wf3mOic7tONWfPl6as+Ks8FFd++hzSdaCh3tm8XkVCbhT4NEAgzLk0bKs6785SOQA0cwbuRmAdJ3oOtF1IpXT6sS77yXDJQwDp5Mi+xyemxKqjbGBUAVPlZ/YWE9F46ag4cRLNRX2C3P2qdHT5dzOS+W8ONpLBcSOE7W5/vL1tdLat9aZo5ZbpUeMF5HBUSUfTU8BJGOFCREWz/ax7H2JN675dCashma+1eQBnRzlqWLmL8WT0/8K4G/61S8AeFNV/RiALzaO8WkAnwaAe7/lTQ7WONQvq3UGGASyDDa+XKbJkGtZY50TXONG7APG4SjxuOTrEQNEvoyhwx8ubp1a5Qjl0UpP8kuYQq6IWOgAwI23uNSN2044LOuwUkbOyhdb4r7RPcoa99tTwHd16wnYAgYafNSfrtsAmrz+KS3yrhNdJ7pOpHJqnbj47W+wzp/SPfwGnT+lJkY+hlDtzCEQK+npt/HJ1zqPqjbIZy3sp0dBl+loal6qmocKQBjTCkDIqVpCpiSnSn7zbfLSLbtJzHVelYT8DI0xzzET7a0yqrGS+pJftYHSEQB7wM3pOfgR2i2AG3qs4j016iCqRLnRb6nsP+m5ashRpIqIPsrMb/nVPwJAenz8KIC/TETfD5eA+AkAP7nooCqMIU9Awh0pmGQgAiSJZqtBhCmEO6DCIUlOSYx3eKSQVRmjx31NecgjBxJpu1zDGMHEPVPGqF5cogRHfGi1pEUACWjocg0W+fptg0f8m88ZqVrimAeMuK1tkZMuWyldJ7pOdJ1I5dQ6QQSYwbrfWyJUsQfbIQzEuTH1PKpcZMJkPbef9lLJQJ/5+cVLVQv9Jcfwf4DkPMXvuUYo5P6P6eBRCwG676CeWyX3lY+2HgigynR31z3iMA4nI1Zyn0eHAwG0QoLu2MuJ1SypIqK/AuBbAXyQiL4A4L8G8K1E9M3++j4P4D8DAGb+DBH9CICfBXAA8N2LenTA349qOAKRDsCgwOQE1jkIcfoOOZkvF6AQq1zckyGvBPBLEUikBatZ6RFYkFgNEa6QvEwz83EslVpSLpCChPsVoMjXkQIHUISCoOodCx5JKMsq8LCxzlSOyTRgpODR3J4/pAnpOtF1outEKk9LJ3TY6DYI1W6I4b6L4VAmpld6++mwn5X8KVDiudJDKEx5qZqhPx82lDK9feErgnv2FF5u24urc6nQ9FQBpbcqD2tuaCwIlTWOSDFbwJpg0D0NYuXOMEOsgOr5l3qtlvT++45K8Q9M1P8+AN83d9xyRyRgAePuSyfnToIIUjY5ByKhDglq6W3pMjvEqIQ/PJAIcDSt9Agc0R0bP2ptSRCddqBD/YFpgEjXK8ABpOBRscQL8MjHQQKOBo8cDJpggbR8iXVe204CMIueadeJrhOyjK4TeEo6QWlCOhGf3EPV7OnXGD5BiISQqVbYrzauVM1LlYuE/mLIL4YPc/KgO1ssFQm7W6+rlggDlIdKhQClCTCEZJT1fIgF8VYlYUCiEA7U50YY4NQTxqdIrICbkauWnMeI6kqZWZcBSe8nYAZEJFwyASIgFeqAWNcKOBjQIQ0BCnnRIfzBenvMLREgKa30GPYAImjoD0CDykkfL+vlEjRkPYCBrAMJcIT9NHhUBi48BXiEHBJGZfBDvx3LAAMeHOYApa1hz0C6TvjlrhNdJ5wQUAyZELxVJyBUO+PIUwj9UZZHNRP2k+T0POynh1Boean0tYXjZiSsNnzC7egGRVKB0lPliFb5HGSIhWQYFMQwIABs/K9lAoYR8OO5+7fmtp2IWIV7maqHmU+c1l/DeZAqIDYK0j08AYrMTl4LIoQwGKIuC5a1USAiqBVazBQoxCon39rluSU5kEiog4N1rsEECjy4cCseEwtvfUQ5YEiZtsBDWbLuHsVkaKOxfGrw0OXAHGCk4AG9D0rgCABzTtJ1outE14lE8kE99fhOawiVhA8LD5UP/YXE9IVhP+2limRq3kuVrDdDf5R4avVxZb0ltST1KdEEKl+W66xJLWk9EEj/nGH8cYzXHWvAIRTIgF8GTkOsgNRAa9bxv82nNOG1qslZkKrQePhfbUWfBESAZJTpUCZHtHDWsm/8Q+hDgEHt4a6xnltSAIm67uBeQHw3FH7TFx9jt8d/TDkYpSEPWcgsdLk39S40cIR9a4BhKexzKvBIwajxq8Em21YDluQPsh/Heusf9a1I14muE10nUpEE9VMQqnw+v3yQzyV5VHNhvzkvlbumeoK6iM6nGtUxynrLn2OtqvQAFH1ME9RdCBDAZMK63E8rDGjYPXceCYOxbpwsY+P1W1MQq+A1u2ViBaRtab3Csus4C1IFIIQzCtCQcqS/KOrOJOrCg4gHCGepKxAJ9XTow++pgCFY6H6bzi2R0EYSIlHrzBTARFvh+j2t6WUwJ8VxEhBJQSOpvwY4VJ0CPDRQHAEeyfpEMi6pfUvA4Xr9CsCEsjORrhP+iF0nuk4AAFQOFcXk7mMJ1W44BEJ14UN/Oo9q8CFAGTE9D3lNhf3CYJ8NLxUA5AnqybE9OavmU6mef3PDKrS8VMwEy+w8x0wJObJI86rkPFPPF2gnreswICxcXlX4dUTK1bUFsQLgZip4isQKmPjkF3j9zoNUSQNQC3OIhY0jQATpQ9KWuAaZGMJgDwaEYKFDdlRnzUEiBxJEIAlAoSx1qBcsgAIAtXDHKSSxyMNCCRqhLsc6VeCQ+tl6beLX3Ho+Cjzyxr1WXgGCVs+mtC6HY56VdJ1wy10nuk54IeBkhEo8U84r5QiTECpDnORR6VHTa2E/PR/fCE142l4qnaAuMmSkbW4sqiWDfh4b/hMPUS0EGBubaW9VLQxoyGLjk9UlvwrWhB6BQqzcjAuu3TGGT0asACwmV8d+/udBqoCo1Dg9iEDqs2vgpcdTsT0JZ2RA4oHCndw3TLXEXQ0k8g2EixZAovATXtxCJn2U1MIc/h5dWQoaoV5SrtYrQJKOCI0UKGplet858BDwgVvOy+es81r9EkjiOc5Guk50neg6EYVwOkJlIqGS35CYnoX8psJ+4qEqxqSywyIvVbjG7B40oZKwXz6yushUzuBaQpUTJEBCfjEEOCX59DVApTcgRtfLkGzIrxqIlYfKulCkcWe31pyMWAFiqC0jVsB6FTgrUlWEOW4KIhT3ARCTYqXHk7eaQ+iDYrl284Xwh5w0f9yV8EY4BlMZBgkXJK0iig9lrtvm5KOsfQWst1NaVljoCgzydX2dDfBIEmmfBngkICXl7Z5NrTor2p+nI10nksfRdSLdDtwtnSCckFB5b5EmVoY4jpQOtZwlj0vYLx+TSienA6h6qQAEL5VcYyv0585nCmKmZcpbtYZQ6Zr5wJ/5uVLP1by3qpbkvhFCNca8KvFQMfkya545sQLS5mmJnA2pIvads28KIgE0fLMfwhn+OBQbSm2JF+16K/yRPN4JK12ti6UOKGDQLzV7t6cKeRQfjT5mDTSkzpHAUeSK6HJp8PX+R4JHfvzyl+sWuzpuHgZx284IQdB1InkWXSe6TgAnI1TS0y8SqzSPKh8+IQ/7jaDJ5PSWl0qGUdDHBOqhP8mnAkpvFfN0aHCNFAOCEid5VXLt+tpAkSTmIt6qPGndgpR3zsJAyBNCz8BAsFQo8HkjVudBqkSRcQIQSepkI+4IwEhbyUgscW2hh2PJ0f2JitySUBOqnrK+laUewSR+EJS/rRUvepFkxyoAw1828vIlwJHVK8Aj7+Gkj3ED8Ji2zjm7BhTXJXXS4yoL/Ryk64Q6VNeJrhPuWZ2SUEli+gAhUimx0mE/IO2JV0tO1+G5OS+VDvnVQn/575S3CkAc9qOxbSphXQhV7plKtoHQugJ3baW3SiR4vTyxAtyky5pMybhhMVEdaBErqrQZx8raPCtgXh3Og1RBg8SRIMKIo0xLHQ8OKLqN++Wahc4pkCThD9mTZ4CkAhxus5wsIlRiCJ7gI6lKDgL5chH6qOwjl77GEq+Ah87vuCl45OesAUPNKs/Xk21nJF0n5Ca6TnSdQPhm1hCq4B0iN7K3S0i3SU+/UA9cDJ+gJfT2U8npOuwnXioAVRKUe6n0bzhHMhJ7OoxC3vNvTvSMBWuk1dNPyFMQ33BMTbbcCgPmwywIwXInMrDEftgFU4QC3VyBdDJiBWC112pKzodUMYC1IEJqnXzDQygTctX8ZREUFEAAykL3QAIV/gjbdF6Iq5sAiSCaXBCndQOYhAtTH6Nsvy1JgKMCKDlohDIPGmo9AYikXrat0mgnDXieZ6LBAfPgkVrnXAWKeC5ZLnNG9L7nhCJdJ9B1outEImsJlR7c05ANPf0kUTxPTM+HTqiF/UadpJ4QqraXSg/2mXuo8tAfgCSnyv2m36FMU8OVbUsl91IlwyowFURoTkauJ7prycOAYZgFIjC758T+WY0ABsMYrZ/7E4A5E2I1JWdDqoLltgZE8nVfqEMaiaWehz6g2nKxkDOwCBa61OO4LZRyBBKwB5PwbjKLN1jmiB+CvhBk9Y6V2sfBleUaaPjfKnD4fdZY4ovAQ9dBCQRLwQNT+7RyRrJzn410nUil60TXiYYISMtI6Tmhynv6SR6VjGIuf3NhP93rT4f8XD0ZQqHtpcrn+Zvq9afzqcK2WT9JlDhLwfKX2OwB6I2f/M70YKC5Z2/KWyWSD7OQe6yQDbUwAn6YhYxYnUhOQazOglSJ1SRAcAyIgBAmnJXGS2/TYYtJC12ABNIAVnJLMpBJLG/yDW9uqcsB87CC1EkeiKp3E6l9bFPAkoOG/HLcrwoc/riT4BEa81qeSQYeK6zxAoRU2EQD01TOiAaim+L2qaTrRP5AVL2bSNeJ51YnWlJMPZNNkNwiVPlEyVNhP0lMz3v7yXY90OdSL9V06C8SKQAJmTpFgromTzmRypPV0/3SEGBOCtd6q3QY0BDD+ndp/btkoCBWMJJbxTDGDQ4azMITeawAHE2uzoJUAYiNAY4DEUABCcV2O2wrljML3e+Qhz9Asb4DDg7Hj/tRABaIRc4oLXV9AeFCZqznUwuXy1XQkN8KcKTryy3xom6WNwKsB4/Q6OfgkSxzA1T09fneUeckXSfSZ3Fb0nXi+dGJTHIPVU6o9FQ2esBN8WYZcJzcOAv7AWkelR6TCkDiqRIvVf0abfYbP7ha6M8d2ygPVRkCnBsYVOTYYRVCGZNzHs0cR4ZXWOqtaoUBhTsl3iq4toKJXGgQCETLWkRideI24liv1dmQKrHOQFhlnQdrXI4j2+Q3B4YJC13vL8sAktySaSs9bk8sdXgwyU6QJPSewAhvSvaxkf5Q5kAjlOXrVIJKXrcGHkX4Qe1/ZL5IFYxmllNQ4eQ+zkW6Ttz8GTal68RzqRNadE5STqiENOU9/WoDfOpR0/NRzvOwn0w/0xpCoeWlqiWot0J/yRhVyDxWTAnQn2pYhVzmxqhKyikdXiH3Vq0JAxrIiO4xv8pVkHdiMFqnr6ZGrMCoGmVHyjHE6mxIlVh0DKQNUgVECtDQ1rgGGKgGHqp+sewbNyrDH+E6gEkrPbxLEtDKLHUuT0z6BI3Gi2ldq0ZTHwBXlsNvBTT8bxU4wvpC4Jho1PPfOfAo9rflugaaslyBhwCJ9UByTtJ1ov5Yuk7cXZ3w0iJUMZ8qHTohn2dP51HpUdO1SNhPj0mlCRUQycYSL1Urj8rtn3m/lEWhc7ZqcpMcIMvAQGhOS1OQKGlQUHrd5Brz59iSVhhwEwYDhfdcuedN5Idd8O6s0ZqzJFbnQaq8AufWOFAHESCCRpp0iyQhF+pP2vG2he6BBBxAbApIala4Pn7YHwgL4b0IoOgL16Le3yQgTEmtLUwAJB63BhqhXAFMup4CR1J/CXjYdJ9J63pheGMSPEIdroKTBqazkK4TqXSd6DrhpUaodG++lEi5bVN5VHNhPz0mFZB6lOa8VLIsx05yqiokTotlSpLUa79zsjZRXY4d8q1Q5lVV9/ENitzfUm9VLQx4sEORX7Xxwyq4RHXG4HsBjklulTQrz5ZYnQepgreILCUgkoc9nDWuQETXIxTj9EAtNy30AlS82U313JJwXiBa6ez2i8ART8LqAsI3pubbKN6TBpdTSHaCRL9y4JDtCjTKshsCR2W/WwEPzsDDAkUSrh5l2p74uZ9Auk6ok3Wd6DqBNqHSg3umRMolpk/lUc2F/ax4qrhMTgemvVR5gjpQjk3lzpcSNX3MEAIEFfuItAC/RqjWEq2pvKr6iOvLvVW56DCgHr/qYE08hxEyhTCGlYUMt+DI1W0RK6D9rEXOhlS5hikFEQYaOSPsJ4DNAEKYqrbGF1jo+jgI2xSQkC6vAIlcf26pQ33Ufp2TE2UNekC5E0uuCzlg6LJayEPfwhRw+P1XgUcOApNgsgI8igRdVsu+zhiXk3s4F+k60XWi60QQGTNch/I0oRKiFIgURbJUTJS8MOwnY1JJeT6EwhIvlU5Qnwr9aWlNdQOU5KolQp7mSJQO/QkJ1GNYxWvSDYMLwdVyrG7irdKTLh8wJPlV8ASKAZcCYBDGsHLPxbU34rW6DWIF1ImqlrMgVb6d9kCwAEQAFwwmgInmrfFamVjQChSK0EYNSDirlwAHMksdEUw8alF+MV44Wz+1FN9BBhhJnRw08rJTAEdlOzGSpN0IDpyt59uXgQdUnXzZbedbfQdrpOtE14muE6XIdCaSNyUTIWvPlCSH68T02nhUOvSnw36tMakAVAf6rF5n5qXS3qlW6K+WTyXnbIX9Fk2v0iABc2SrNbJ6sz4ohDrnvFVzYUAAYZgFeI/VASZNXNejrPtegvCE6vaJVXvbWZCq0CgEIGiACKEyxxmDLS2zxhtlYNWe5wCBCpDk9RAuprDUayGOAiw0sGg59jto6UENLHT9Gmj43xxcWsARjr0CPHT+xinBAx4U0rr+eKNa1iBzLtJ1outE14lCaoRKlmXoBOnppxPTJY9qavgEAEkvv6nkdD0dzRIvFVAfmyoeryRSY1a+JNynZW0uFVCqifTum8qrykOAS7xVcxIGBYUjU5JfJWRKEtcNACZKvFjmKRGrlpwHqQKaIALSZUjDGaquzinRANMskwIFAsFCl+uZARJ5TzH5Fqmlnt5ebEj1Nm+5J3UTFDpeqjpVAw4sAI2i7ObAMVmmG/sjwaPs5cTFsgYZGs8MQbpOxLpdJ+68ThAQSFQ6n594ptKefltFrCSPKh8+QUst7JcPyBmGO1jopQLK4RO0l0qH/mr5VCI6DHgbwyjkg4HmHirJqwpCQC0EqL1Vc7IkDAjEYRYkwSpPXJe8KhnDasSzJVZnQ6p8e1yASOgyXgtnSKKuNOjkWOtia1wdT7fZCZDo9QxIQtiZyjBIeH85aKgyDSphmxz/VFI5FtW268a+Va5zS6aAo9hvGXik1vktgock4Op6Yv2veLS3LV0n4kV0neg6AUKcFFkN7ilDJ+gEdSEySQ5VJY9qadgv1lnmpQLSYRTkXNWcKpgqUYp5VXGb9lBxqLf2Qc5La7qauX1qCetrvFWtMODGWBysCflVOnFe9wgcLZ45saoPfKGEiN4kor9LRJ8los8Q0Z/05a8T0Y8R0S/439fUPt9LRJ8jop8joj+06Epyl7RKntRluuGoNUb5XFbNxkvOYfPzVepwWQ/qWsiS++N4TGTXq/dDdi6Pfen6qf4qx8WSa2RM359u9BvvRb+z5FyNsgI8Ku/mtsAj1FsgXSe6TnSdePo6QeBAqPTgnjpBXSem6zyq1mTJIoFIoT4mVZ6cnofmRMRLpb01rQR1dzxKlkPie4XK3oaHihuErbbeklb4conURplPtzuSLISVKHon9d9grCdWrg4B7pfgh1zgsI6VZHGtzJIqAAcAf5qZvwHAvwHgu4noGwF8D4AfZ+ZPAPhxvw6/7ZMAvgnAtwH4n4homDyDUm6qNUAVC63aCOnjWJ5utHJwmQOSvJ5aj9fkTPeksVUNdnmdFWC5jb/8uVXAMX1WlP2V9z4HHEst8bKMm0QhL78ReAQwlS7lsm2RwnWd6DrRdSKV29cJIBCqwjNFQqzGJI9KD5+wJOyXj0kFoMip0pMmt7xUQDuHKpxXjXuVj1GlpZW8vkRuOjnw1HQ4S3ontgYsrXnB8nG9al6/jSdPQqT0NERCumSQUAJgjA3ESv/eJrGaDf8x81sA3vLL7xLRZwG8AeDbAXyrr/ZDAP4egD/jy3+Yma8A/BIRfQ7AtwD4ianzkNdfsgwY13k2hCGsr+CC6uk2/wtfz+fV+W31vBIwmt3KwbFcjgXZV51wNm9ET8GhbxBIvY/UeL830YWp70VtS86rwxi6Hqu6qiyvp7cfCxzNbRXAz8nGavCwHMEjJP4uU7SuE/45oOtE1wl/y09BJwicEKq8p1/wXlHsEdgaPqEV9svHpEo9SeW0MbmIl6qWoD434KfkU+kk9XpYkKrLNQLV7PlXLUUz7GdBMIyEIrbGrVqSsK5Fcqvq2+I+hl17l+dXGblu1SNwtC7XSkZdt+pXt4S3EQpclVNFRB8H8LsB/CMAH/GKBGZ+i4g+7Ku9AeAfqt2+4MvyY30KwKcA4OLeq06BoUCEADZUAIXbRknDD/+c2KhGPvz5Y6m8Eg0yQLnPHJDofSbzRkgd27+8kHMiB6l8Zzrv5CZS/Yazj6gADLWcg0ZRlq3XQKJVPgseNdBogIezzE8AHm2jsildJ7pOdJ1I5bZ04v5HXspCfdxMTNd5VK1R0wFPZJAmp+dhv6VeqhaRmpOlkyPXPFanEJmqRuQYNSt6AGJ+eIVablWetK7L9WjrAMAq10r20z0CrX8vlumpEqvFpIqIXgLwVwH8KWZ+J0yGWqlaKSveEzN/GsCnAeDhqx9jslwABhRYpGPy+GRdWRdr3MbnUwAPlIWeA4Y2oNW2GlgEIBHgUCcJx1HWeqjr1yN4oHiZSdLvCaU4XgUskno10KjVrQCEXr+RJV75LbqS18BDH3cGPJLr0HWXPteuE10nuk6kz/UWdeL1b/gQ66ETaonpOo+qNmq6Fgn75cnpuUhy+pyXCiiHUdCSe6lyMqXzqXSSuiTE167rtmTNGFU5aWrtu9ZbVZ0bkFxvQADOS8WeyImXCkh6BOZDMAB068RqEakioi2covwlZv5rvvhLRPRRb318FMCXffkXALypdv8YgC/OnsMCTpsBkB6Px0OFRdLoi7Wd1s0ae6rvw17RC8Dwv1WAYXU8qO0ogSa8I1bHg9oXlXeo692GZMdOztUCiImyZDtXts2WLwOPmB+S1Q9W+2mscbe+/AV0nUDXiays68Tt6gQhDp2QJ6bX8qi25rAo7Ncak6o2hMKUlyonDHqcrJxkyFAN+lrWyJrE8GPm/5NjtkiQHl1dJ+WXPQCFFM17q/IhFrS0Rlt3QUB3PhkYlOGbRmJHwBThsv6ib5NYLen9RwB+AMBnmfn71aYfBfCdfvk7AfwNVf5JIrogoq8H8AkAPzl7JVljkSZgVnIFuFK34h4vGh3fWEjSrm6okkTVrCw0Ztk+SWOGxr5Z8mvr2Lf6V7vH/Ho42y73md+r2j955rLN1p+7NNRkOTlO3FZ/bqcGjwTENHgsbHe6TnSd6DqRytPQCUL0BOWJ6TqPSg+fsCbs10pOXzpGVPSMTSeoa0mGbFD5VFOyJvl8ad1WPVa5XRYU1mtSEL2KM7J2b0uS1nV5HEHfhuVa4voQktoRkttjMjtuLXl9iafqDwD4YwD+GRH9tC/7swD+WwA/QkTfBeCXAfxRAGDmzxDRjwD4WbgeId/NzOPkGUTBCQDEylbj8WRW+FToQyzh3KrWIz1HC54RrHRC3Ur3pw7H4Wwb6QP6+xFrXcrlGFB1UdZJ5CakeeLbqFrjrTr62ifKdKOs18vyDAhqIF8hCClopOBRgMwMeDStcXWMBdJ1outE14lUbl8ngNDTTyemB48QlcMnrAn71ZLTgeO8VFpqCer51DOtqWniMAvThKYlMvffEjlmbKq4bzsEuCaUuCYMCO/VCp4oQpK4bojB5LxZMoaV+KaAdPwq3Wre1GO1pPffP0C7OfuDjX2+D8D3rboSaWSsuykeUmCIDf906IOpksgrh2cU4Y/FQKL2ZUQgCe+AyrIENFAHi+Jbkyd93LddysxxaoCSA0itPPdChO0aGELdGeCQeglI1PcJ4FGrewNr3E0mqy984pF2neg60XUifaRPSSdi7zou8qj0NDRLwn5TY1Kt9VIBsdu/XFfLY5WPoh6PXx+fak6mSNPUtpuQqOQ4oGYIsJawXsutmkpar5Xp0daBNL/KqnqSW2WBkNcpg4MCt0Oszm5EdcDdFkbfkhuGtsbdKNIORHKQCC27rl8BDDnfWiAJwKCsb22pA0itdagLRKxTLEPteyrg0KeqHZMb2ydAI6w3gGQOOJL9m+AwDR5lYq5fV/OWtcDD9ZCaAI8V+SNPQ7pOdJ3oOhGFKE6crAf41HlUa8J+tTGpYr2YLL7ES7V0ahaRnEyl5z6OXNVkjacKQJiOZgnZynv8SdncUAtzor1VU70BgXSYBQsVLqzMEShDLQAAKQvv1MTqPEiVKLtv+PNQx5LQR2jMoR/LdPgDWAgknB4/PG9G0W1cW+aAei+5la6+r+RabknmQGQSMHTdCpBMAkcowzRwqH1j+EL/1nKIGDSm54xAgfXgcUYA0nVCXcstSdeJ50snCCgS0/M8qjVhv9qYVLWBOFu972qSe6laoT/5TYdqoKLe1LnZhwRbIgnqx06sXCNCeryqYtuCXoBrvVVaZsOAANi3UwCKxHU91IKMYSWKFjup3pxYnQepApyyG9/sWz+HmIDEktBHBSQEZPLwh667GEhUnRwoqpY4JgBFf0saWJ6GqHNVrXFdPgEmSZ2wfRo4wjGqIMKNchTgESzzLHekBh5z+SIJeNin/C7mpOvE05GuE8+NTrTyqNaG/ebGpNJeKpEpL1U+z19L9Cjq+XHd9jqI69wqy4Qlr2QJmZrrGei8RAi+vuZ0O1gWAlwqNW9Vur0eBqzlV8VrROgRqMewArHyVgE3JVZnQaoIrjFgKBBB/HAmQx8eKJzFvCz8oa3rOSDR+zM5ZGmGOuT5+/XECkdjufKt3TBPLt5HTbLyKog0wKQGLGHE5RsCR7VObYRolSuSAIMGmmPBw6YW5bOUrhOpdJ1A1wlw0ttP51HJ8AlbirnuU2G/1phUuazxUi0VPZSCPk9cXp+Q/iylNYxCa3mpt2oqaT2s+/0lDAhU8qv0dcgo66AkRHhKYnUWpAoQJSenzAT3JAhg4y1iNEIfiG18DhQaLOTxuIXlQJLWYQUUE2Di1wGkVni+PWvMbyvkUT1eDTh0Xa6vx/Ky4U+WG6CQlnGybcoSr4JHkkuiwEN3JV8DHmcU6gC6ToS6QNeJrhMA4kjlwTuFMo/qJmG/m3ip1ob+3K8pzqHrnnoEddb6GcroaAUrCVV9PKpjvFVatLdqKgwIKvOrdOK67hFoQGCKvQVBDDcySN4grSNW50GqRLE1KPjQB/wI0nnoo2ahCw5MhT/crguBRD7A3PoGZsEkeS95GSrv5/jveplkxy7OlQOG3qcFGn5bAjI3BQ6gAI/J0IYCDwcKNwAPi/MBka4TXSe6TiRCwOI8qrmw35xIcvopPUZ56K+WT1VLUj/2GiS0d+zgn3PHzvOqpkKAedkx3qolYcC5/CqduJ4PtaCns7kJsToPUgWv/N4Sd+BAAURqoY/EQpeHNxP+yAFiEZD49VbyrVsuwcRfaAQdX1a1vCsW+q3IFIhUlnUYI6mfAYKUFaDSAI60TP8uDG3k+4xc1mcGQvky8KAzAQ+RrhOnfZ5V6Trx3OiEDODYyqMSqYX9ZLk2JlXNS6VlqZdqieShv5YnKicPx45RpX+nxLLz3MxZMdY3CLnXaWkI8FTSCgPWhlkApflVjJi4rodaEGJFeQMFCh6+JcTqbEgVfKiDjICD3IVvnAO41C30PPwRQYDT/ZJtE0Ci6mjgAFIwSaxuIH6UCWjEl1BMhfW0wEMkP5dazwHDlaX1pizw1vIa4KjttyRXJKnvQSIBj9Gq9QZ4nFH+CICuE09Luk48NzqxNuwHRDI1NyYVgGQIhbUk5tjQn/zmPQBrU9fMkaubcmBJTM/JVc0zleyTEKoYAqwlrJ/KW1ULA8rwCyAVBuQ0vypcJ9IegVDb3dXHhsjlW1EkVhNyHqSKEbqPMxSIoAx9FBZ6sIg9OvjwR9m76QggQQSSBFQETCpl4dvw332w2gHkYJFY8FrW6XJbGu++sD6XAkajvL08AxwAjrbEc/BQ3cOPBo+Vk8feqnSdSKXrxJ3XCcL6sJ+eiqY1JpVIjcTkXipg2RhOucwlxS/1WM2Vn0o0SQEmev0t9ELdhrdKy1QYsMivUteT9wgE0lHXa8RqTs6DVAH+ul1CZUim06GPCQsdAObCHxE0FgBJHjoJKKTrxLLEWvd1qoCSgEelIZf9T/3tVY6XfN81ANHlDTBpW+Nt4IhlDeDwxz82V0SDRzrQoQOPsru4Bg91kecgXSe6TnSdSGTJ8AlADPsdk5w+5aUKXjG0E9RbknjEVD6VvmbZlnqsng2xmpNarlY6gnpjnCtwse2m3iqRWhjQXSuHtiRPXNc9AoVYUaOhWkKszoNUsbOceHAcUezxJPRRsdALMGiEPxJLeymQIAIJEyXWdwA1DRLaWpeGUj97tZy8L739NtuuFSBSAEZeVgGU2W7kYd8MOPT+C4CjGtrQ4JF0GVfgMXqrvWqZS9kZhTq6TnSd6DqRiLMd6qOma9FhP5mKRhOqpaK9VGtkKvSXXB+o2Nb0WDU9WdPXchtJ6nI9U3lVrRBgcZxGb8GarEtarw+zEI0+uc6yR2A+hpWT2BgV6QqZnAepglf00bpkXELsYcRoWuhMWBT+WAUkGiACYDAiQCgwwcyvtKkaHHLwuC1LvCX5edR6EzAm17kNLg3g0Ntz4EiOuQQ8arkiR4IHnVGoo+sEnt676DrxnOgEVof9alPRHOOlukmCen5u91vmU+Vyql6Ap5apXn1r99NyzCjrQOqtSkOWlWEW2D3XMHUNYnK6JK7XxrBy4hotmmFVZ0Kq4BoBSO8mgOUGGhZ6AIOZ8IdriFYCCSKQ5GGNJpgwUDx/VMrUN8PZelXW6tHM8Ypvmhvbc8BIyrjczq3lBcCBtGy1JZ6BRzWcUQ115OBxRugBdJ1oSdeJO6kThNOF/eZkjZcqT1BvHjML/S3Zntet7Ts1Vc2pvVRLCNXaEOCct2ptGDA933x+FeBUTPcIrI1h5cQNu1A2GlHOglQ5i5wB41mgtOATFrpLup0Pf4COBxIIOOQAUQETtx6Qqw4oGXhU30uuHzfVicb+NWBola/pRh4a82p5Bj4TYY2wzxLwCAm5aa5IAh6FZV4BD1YX/Iyl64R+GNl614k7qRMgrAr76aloanVv4qWaklrorzaKuhCFWmhPJ80vmRC59YZOFf6TXoH5sfLjHxMCvC2J15ESrTy/Cuzay2QEdsSmq0aspt7IWZAqwCk6W3L+OWDaQh8jkMSwhAOPJPzhAWI1kHgAi8sKGBpg4tYbgIJYL9atPQTcTvtVOeaUdb64G3kBKDxhpU8DRzzmcZZ4FTzy3kwaPPy2BDxGezb4AXSdCHW6TnSdADzwLQv76TGp4vZyCIWaLPFS5QnqLamF/vQ5ynDg8QSqqNcgVGuJ1ikJUex1t3x4BeA4b1UtDAjM51e1hlqwzUYqlbMhVbDWW+S+wZ+00JeGP1YAibbGpSEPgJGCSQIeGkyAElAQ67lyKt+L6NFtN15VIOHq9ipgVNc5rlfBZQI4IOsogQNp3UXgYW0FLLL10YZlUpY5zizUAaDrBGKdW5OuE8+VTrTCfgCSsF8+JlVNjvFSta9rWehPzqt/W3WX9ACckinitIRYCclonXnpJMq3MvjnkWHAMK4VUB0YVK51DHlW64nVeZAqaRyCRY6mhR4abN+4N8Mf0sgcCSTV5FtSxwRSaz1ceAVQ9DvIX8htWeItyU/P9W31kIesc7ZescBluZYfovdRwBGOuRY4WmBR6x4+BR7nFOroOvH0pOvE86ETjeuoTUWzJDm9JnNeqrkE9Tz0J9enr6PYPhkGjGXcWD61tEJ9RR21fU2ieqvuGm/VGmmFASX0lw8MWktcL4jVhJwHqQK8Wzq3yEsL3YUygCXhj8VAEkIUFMFA5lDjFExSizwLh2jQQLlckNwJ8Lipzkx+32tBpAIY9fUIENWeS75uAhzAdFgDqu5S8GjlimjwyHs6MQPjeD74ga4TuXSduNs64T7XdtgPiN6nJcnpa7xUNVkS+qsNpeC2pUSqGG296R9aL80wICKBcud0npmlHqwl26eS0FtjVi2VNWHAfNJlkWL8KilX6zViNSVnQ6pCIwVkFrla54XhD/B6IEFmkY8oLHUIUMiFJb/uelOLPILKEvDQen0qb2nzOJMgwmX5FGj49cICl3p63xngCMe8qSVeAQ9SYY5k2zg6MqDu+yyk60TXia4ThbTCfnNjUk3N7wdMe6n0MAq165kjGVP5VHPSIlhz3iohU6fsAVgjaLpsLgQ4T8hO461qESsJA4KUR4qREKtW4noyOOjEJZwHqWJ2cX0iwFDDIicwaWDx4Q/kQOIb4oVAEpYFSBhqwETEVr0CJsUAiL5e+OUKqACpBS+PIK9zYqkeO2swp8Mbup5q7G1eP25LQyATIQ2o5QpwhOOuBY9arogGD2uDdc5Sdi5medeJrhNdJyZFvFRueX5MqmL/I71ULRI1F/pz11OZCidLUp9KWn+WY1VNDY9w7LFuy1tVryvHt7A8pMMs+DYr9bJRlViNM5d5HqQKiI2Z8sEVFjpzCSyU9YhiXgUkLnSSAomEVNIwCBIwARCuS9ZrgFK3xmNB9Ajc8PktkDqINLbngKHKqqAh6xp4xNOyAjiSfcO24y3xAjy0ZZ6Dhz2fUAeArhNdJ7pOZFIL++kxqfIJk7XkQygk2xZ4qYpyaudW1UJ/tWT1pTlCxyaqn3RIhZxQcTuvamkIcEqO8VYtCQO6kCfHNgvT+VW6bAzPs/1Mz4NUMdyUHMab1mJiFxa6qx6ABVzvEbUQSELoQoAE8Mm9CGEQoA4mOhwCIFrsQPKyii7k6l3ctiVek3kQ4aI80Yc50ABQDWf4ukcDhz/uavBIAEPXU+DhQx0CcmchXSeemnSdeD50gpBeiw77yZhUezskdaaGUJjzUtXKayRqSYJ2ct2VvKlWblUtxPe0vVWMVH2PlVYIcO1goLlob9V0vXoYMFybIlo1YqV7BE6eZ+5CiOhNIvq7RPRZIvoMEf1JX/7fENGvEtFP+78/rPb5XiL6HBH9HBH9odm7BUrrSSu8tSBr3ci/o1rXDZDfn0b2M7C7hgjWl4djwdVR82G57a4Bc/u7hkv2Cdutf9p+W9yXQSPcn0VYDvXyunJu2/jjE/1Vji0Nvr4WfW0OzNN7Kepb/9z19epnNfrGf+TYZTy5b47vktWzDs/Slpa4vBcJTxTr/s+v02gnwYOtBY82BQ9WmcRdJ7pOdJ04P51APezXGpNKSz6EQihvgGRrsM9agnor9Jd7p+Z6+E3J2tHUAZw0n2rqWji5z5IsLjlGTVoeqVbnASAluVXPo39/xk/OLWUU3jfHdfWuSZVPXfkST9UBwJ9m5n9KRA8B/BMi+jG/7S8w85/TlYnoGwF8EsA3Afg6AH+HiP4lZh7bp3BKDSKQt8hLC90fHzWLPLPafSM5aaWrnlEuRMF+m8oxkRMaZbkzkFjrvk5ukQerXZWFu9UWfPLw4uM4mTSOVbQBql5rnJ6a9R3WdSjD1y0scCAAfnKMZPtEWEPqrwltyPZargjbuH0cXfmyZ991outEstx14mnoxLKwX1m/DoF5eQTUZV6qpaG/2rmm8qVaXqtnIbWwXzGK+oK8KvE+Lc2duqm3amkYUMrivSBJXBdSR7IvTtD7j5nfAvCWX36XiD4L4I2JXb4dwA8z8xWAXyKizwH4FgA/0T4JnBKrq5Xej/4iUIQ/iCZyS/wxJoGk0u1cHqYOZ/jGLIZCMjABUkDxx0gARbfHGljCzU48zVNJC0hmwMLVkbIMJCqgEY6pl23cfw1wxGNVwGQpcChLvAoe4whYKV+mxF0n0HWi60T6KJ+GTnjJw36t5HQtc16qnBS0BvvMCcESgjCVrL5mJPWnLcfkY7Xyqlr1liSsL82tOkUYEEDIrwK8dyq/ZmCyfVqVU0VEHwfwuwH8IwB/AMCfIKI/DuCn4KyUr8Ep0j9Uu30B08oFiFWuQKKw0EWEKnLLIp8HEig8Cvkj7JYlx0Rb6iHPBIjWoVzSBKAAQKIrlDXYcpx13+3NpHauGlgAJWAA06Dh13MLPNSv1smAw+9XBQ6oeseCRx7a0OChr2GhdJ1A14muE4nclk4wIiHSU9HkY1KJTA30uTSXqjaMwpLQn5wjD4VpElALbS0hBXLsJXJsonp7gM4jBv1c4MmqX8NxCe4iLW+V3i5XqIdZkPpT+VWT5116gUT0EoC/CuBPMfM7AP5nAL8DwDfDWSh/XqpWdi+ugog+RUQ/RUQ/dW0fpwqs8kaSXACdK9DKLQkNSq1M5TZIo5PkiPi6kpug80z0Powy3yTLO6HsuLV8DueNeMp/letI81qya9a5ICNX90n+RtfgY9T3Lc8ZaX6IPF89H1mWI5KAyVSeSC2PZBxdCKOVK2IVqIwjeBxXqX7Xia4TXSeenk6899X9oqloamHA2kCfc7lUSdnK0N/cUAq1JPW0fnqNkrM09S5q+VW3lVM1d/5Wx4DWNiDe69So87nkHqwpEqbfqRCjPD8uyadq5Fc5x337uS7yVBHRFk5R/hIz/zUAYOYvqe3/K4C/6Ve/AOBNtfvHAHwxPyYzfxrApwHgle2HmSXJlggw0t2I6uEPqZfnllQt8qwsWPN+d7Af0ycNg4ilHUIhiMcMOiCGuu427gexCKERII79o+9fW/Hpw64UHikTFmZxbsv17azLOClLrG91vlkLXOq2rPBamVXb1lriea4Ie0AJwGVjvYUQ0nXCX0DXiaSs68Tt6cSb/+rLvCbsNzUdjRYBz6K8kqA+tV6T1tAOc8unEPFSrfVWterfZGqaqW3Tdea9VUuHWJjaVhtt3QLFwKAnmaaGiAjADwD4LDN/vyr/qI+jA8AfAfAzfvlHAfxlIvp+uATETwD4yemzOOUN9yRx/ARI8vCH7MppHcwBCSEZiZrINfpZGMTdIwow0eGQmGMiLWry4NTtVcCj9hgqQHMSsfVjFt9xDSxUeVKmQCPPL0lzR2qgwmH/RcAh+5wAOIrQhgIPXjglR9eJrhNdJ7Ln9lR0Qk2avCI5fcpLlXovysE+8wT1fH1p6A9Amc+VEcJ8een8fnMgfwpvleIW9WtQIb5TDBK6Jrcql6VJ63q7XGWLWOl95ojkEk/VHwDwxwD8MyL6aV/2ZwF8BxF9M9zz/jyA/wwAmPkzRPQjAH4WrkfId8/16HAu/9E38CYm3lrrQCTcHSAWeZFbcgyQgN0TC8m30rhnYCIWdA4myCx2aetyUFHl7hrUvev3OzdU6ymkcopWYm4LMNy2dL201GvAsgI4/DGqwCHHOgV4jBZgC5Y6y6TrRNeJrhOp3L5OqDyqNcnpc1J4oFAfRmFqXV9HbfT2/DpanpXJUGBxf9VDBDkulyo0FZOyxPtVHzA07QXYSlhfK2smW25PYQMAbrT1Yh9Jvsd8zhSdw9xORPTrAN4H8JVnfS3PSD6Iu3vvwHnd/29j5g8964voOnFW38SzkHO6/64T5yHn9E08Czmn+2/qxFmQKgAgop9i5t/3rK/jWchdvneg339L7vJzucv3DvT7b8ldfi53+d6B5+f+F/f+69KlS5cuXbp06dKWTqq6dOnSpUuXLl1OIOdEqj79rC/gGcpdvneg339L7vJzucv3DvT7b8ldfi53+d6B5+T+zyanqkuXLl26dOnS5XmWc/JUdenSpUuXLl26PLfyzEkVEX0bEf0cEX2OiL7nWV/PbQgR/UUi+jIR/Ywqe52IfoyIfsH/vqa2fa9/Hj9HRH/o2Vz1aYSI3iSiv0tEnyWizxDRn/Tld+L+j5GuEy/2N9F1Yr10nXixv4kXSifYDwz3LP4ADAD+BYDfDmAH4P8B8I3P8ppu6T7/bQC/B8DPqLL/HsD3+OXvAfDf+eVv9M/hAsDX++czPOt7uMG9fxTA7/HLDwH8vL/HO3H/RzyvrhMv+DfRdWL18+o68YJ/Ey+STjxrT9W3APgcM/8iM18D+GEA3/6Mr+nkwsx/H8BXs+JvB/BDfvmHAPzHqvyHmfmKmX8JwOfgntNzKcz8FjP/U7/8LoDPws1Gfyfu/wjpOuHkhf0muk6slq4TTl7Yb+JF0olnTareAPArav0LvuwuyEfYz4nlfz/sy1/YZ0JEHwfwuwH8I9zB+18od/n+79w30XVikdzl+79z38TzrhPPmlTVJuu5690RX8hnQkQvwc1g/6eY+Z2pqpWy5/7+V8hdv/+avJDPpOvEYrnr91+TF/KZvAg68axJ1RcAvKnWPwbgi8/oWp62fImIPgq4mdwBfNmXv3DPhIi2cIryl5j5r/niO3P/K+Uu3/+d+Sa6TqySu3z/d+abeFF04lmTqn8M4BNE9PVEtAPwSQA/+oyv6WnJjwL4Tr/8nQD+hir/JBFdENHXA/gEgJ98Btd3EiEiAvADAD7LzN+vNt2J+z9Cuk44eWG/ia4Tq6XrhJMX9pt4oXTiWWfKA/jDcJn+/wLAf/Wsr+eW7vGvAHgLwB6OYX8XgA8A+HEAv+B/X1f1/yv/PH4OwH/wrK//hvf+b8G5Zf9fAD/t//7wXbn/I59Z14kX+JvoOnHUM+s68QJ/Ey+STvQR1bt06dKlS5cuXU4gzzr816VLly5dunTp8kJIJ1VdunTp0qVLly4nkE6qunTp0qVLly5dTiCdVHXp0qVLly5dupxAOqnq0qVLly5dunQ5gXRS1aVLly5dunTpcgLppKpLly5dunTp0uUE0klVly5dunTp0qXLCeT/BxEUvdH6Hjm9AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "edt_1 = distance_transform(image_1)\n", + "edt_2 = distance_transform(image_2)\n", + "\n", + "print(np.min(edt_1), np.mean(edt_1), np.max(edt_1))\n", + "print(np.min(edt_2), np.mean(edt_2), np.max(edt_2))\n", + "\n", + "fig, ax = plt.subplots(ncols=3, figsize=(10, 3))\n", + "ax[0].imshow(edt_1)\n", + "ax[1].imshow(edt_1)\n", + "ax[2].imshow(edt_1 - edt_2)" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "ac57d32a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(2, 256, 256)\n", + "(array([0, 1]), array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,\n", + " 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,\n", + " 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,\n", + " 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", + " 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,\n", + " 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,\n", + " 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,\n", + " 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,\n", + " 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,\n", + " 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,\n", + " 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,\n", + " 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,\n", + " 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,\n", + " 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,\n", + " 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194,\n", + " 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,\n", + " 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,\n", + " 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,\n", + " 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246,\n", + " 247, 248, 249, 250, 251, 252, 253, 254, 255]), array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,\n", + " 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,\n", + " 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,\n", + " 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", + " 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,\n", + " 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,\n", + " 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,\n", + " 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,\n", + " 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,\n", + " 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,\n", + " 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,\n", + " 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,\n", + " 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,\n", + " 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,\n", + " 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194,\n", + " 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,\n", + " 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,\n", + " 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,\n", + " 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246,\n", + " 247, 248, 249, 250, 251, 252, 253, 254, 255]))\n" + ] + } + ], + "source": [ + "dim_for_interpolation = 0\n", + "values = np.stack([edt_1, edt_2], axis=dim_for_interpolation)\n", + "print(values.shape)\n", + "points = tuple([np.arange(i) for i in values.shape])\n", + "print(points)" + ] + }, + { + "cell_type": "code", + "execution_count": 189, + "id": "cbb8630d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(65536, 3)\n", + "[[ 0.5 0. 0. ]\n", + " [ 0.5 0. 1. ]\n", + " [ 0.5 0. 2. ]\n", + " ...\n", + " [ 0.5 255. 253. ]\n", + " [ 0.5 255. 254. ]\n", + " [ 0.5 255. 255. ]]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "def xi_coords(shape, percent=0.5):\n", + " slices = [slice(0, i) for i in shape]\n", + " xi = np.moveaxis(np.mgrid[slices], 0, -1).reshape(np.prod(shape), len(shape))\n", + " xi = xi = np.c_[np.full((np.prod(shape)), percent), xi]\n", + " return xi\n", + "\n", + "xi = xi_coords(image_1.shape, percent=0.5)\n", + "print(xi.shape)\n", + "print(xi)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27bf652a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 236, + "id": "df2d90d0", + "metadata": {}, + "outputs": [], + "source": [ + "def slice_iterator(slice_index_1, slice_index_2):\n", + " intermediate_slices = np.arange(slice_index_1 + 1, slice_index_2)\n", + " n_slices = slice_index_2 - slice_index_1 + 1 # inclusive\n", + " stepsize = 1 / n_slices\n", + " intermediate_percentages = np.arange(0 + stepsize, 1, stepsize)\n", + " return zip(intermediate_slices, intermediate_percentages)\n", + "\n", + "def interpolated_slice(percent, points, values, **kwargs):\n", + " xi = xi_coords(img_shape, percent=0.5)\n", + " interpolated_img = interpn(points, values, xi, **kwargs)\n", + " interpolated_img = np.reshape(interpolated, shape) > 0\n", + " return interpolated_img\n", + " \n", + "def interpolate_between_slices(image_1, image_2, slice_index_1, slice_index_2, interpolation_dimension=0, **kwargs):\n", + " if slice_index_1 > slice_index_2:\n", + " image_1, image_2 = image_2, image_1\n", + " slice_index_1, slice_index_2 = slice_index_2, slice_index_1\n", + " ####\n", + " # possible extension, handle all label ids separately\n", + " label_id = 1\n", + " image_1 = image_1.astype(bool)\n", + " image_2 = image_2.astype(bool)\n", + " ####\n", + " edt_1 = distance_transform(image_1)\n", + " edt_2 = distance_transform(image_2)\n", + " values = np.stack([edt_1, edt_2], axis=interpolation_dimension)\n", + " points = tuple([np.arange(i) for i in values.shape])\n", + "\n", + " output = []\n", + " for slice_number, percentage in slice_iterator(slice_index_1, slice_index_2):\n", + " interpolated_img = interpolated_slice(percent, points, values, **kwargs)\n", + " output.append(interpolated_img)\n", + " output = np.array(output)\n", + " return output" + ] + }, + { + "cell_type": "code", + "execution_count": 234, + "id": "17bbcb9e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(7, 256, 256)" + ] + }, + "execution_count": 234, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result = interpolate_between_slices(image_1, image_2, 30, 38)\n", + "result.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 235, + "id": "ccdbcd2a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Dims(ndim=3, ndisplay=2, last_used=0, range=((0.0, 60.0, 1.0), (0.0, 256.0, 1.0), (0.0, 256.0, 1.0)), current_step=(0, 128, 128), order=(0, 1, 2), axis_labels=('0', '1', '2'))" + ] + }, + "execution_count": 235, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "viewer.dims" + ] + }, + { + "cell_type": "code", + "execution_count": 191, + "id": "5e831a75", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(256, 256)\n", + "[[False False False ... False False False]\n", + " [False False False ... False False False]\n", + " [False False False ... False False False]\n", + " ...\n", + " [False False False ... False False False]\n", + " [False False False ... False False False]\n", + " [False False False ... False False False]]\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 191, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAD8CAYAAACVSwr3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQTklEQVR4nO3dfZBddX3H8fc3m2QzCaBJAzGEKBFCa1ABGwOCWloGiVQJ/qETH1PLTKgNrQo+BGxH+wcddESkU1GDMFJrCQzokFZGBYpjLSAmlKcQQyIEWRISENTwYEI23/6xJ7rkt5u92b13z73L+zWzc+/97e+c88kBPrnn3HsOkZlIUn/j6g4gqf1YDJIKFoOkgsUgqWAxSCpYDJIKLSuGiFgYEesjYmNELG/VdiQ1X7TiewwR0QU8CJwK9AA/A96bmQ80fWOSmq5V7xgWABsz86HM3AmsBBa1aFuSmmx8i9Y7C3i03+se4PjBJk+M7pzElBZFkQSwnaefzMyDG5nbqmKIAcZedMwSEUuBpQCTmMzxcUqLokgCuDmve6TRua06lOgBZvd7fRiwuf+EzFyRmfMzc/4EulsUQ9JwtKoYfgbMjYg5ETERWAysatG2JDVZSw4lMnNXRJwD/ADoAq7MzLWt2Jak5mvVOQYy80bgxlatX1Lr+M1HSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVLBYpBUsBgkFSwGSQWLQVJh/EgWjohNwHagF9iVmfMjYhpwDXA4sAl4T2Y+PbKYkkbTiIqh8ueZ+WS/18uBWzLzoohYXr3+dBO2owblicewc+rEIedNfvg39D7w4CgkUqdpRjHsbRFwcvX8KuBHWAyjZsfpb+R9F3+PpS/bPOTcMzYs5JHvn/j715OeSqZ//fZWxlOHGGkxJPDDiEjg65m5ApiRmVsAMnNLRBwy0IIRsRRYCjCJySOMoa55R/HgZ6Zw9jH/3VApAKya+32Y+4fX9+78He8/8FwO/eJtLUqpTjHSYjgpMzdX//HfFBE/b3TBqkRWABwU03KEOV7SDr7t5cw74B5unL5+ROt5/cRJfHPZlzlr58eY8S+Ww0vZiD6VyMzN1eM24LvAAmBrRMwEqB63jTSkBhbjx/PaNeP498N/xAUjLIU9/rR7Iv/1iS/wTw+tYdyx85qyTnWeYRdDREyJiAP3PAfeBtwPrAKWVNOWADeMNKRKXTMO4YBbX87FM+9q+roPG38AJ0zqYtX3vsW41/9J09ev9jeSdwwzgJ9ExD3AncD3MvP7wEXAqRGxATi1eq0mGn/4K3n88qlcd8TNLd3OhOjiX//zG+w8bX5Lt6P2E5n1H94fFNPy+Dil7hgdYfysQ3n0Ky/j3gVXj9o2Vz07mYv+8UMcuPKOUdummu/mvG5NZjbU8n7zscO88MqDR7UUAM6Y8hy/etdzo7pN1ctiUEOuWXA5T3zkTXXH0CixGDpI10EH8Y1rvlLLto/t7uaHF3yR371jQS3b1+iyGDrJuOCV4w+obfPTu6bw6yPHExOG/rq1OpvFoP1yz6cu47GP+inFWGcxaL+t/vilMK6r7hhqIYtB+607JvDISr8VOZZZDBqWu068gp7rj647hlrEYtCwTB43kWNe0dhVnOo8FoOGbfbkp+k66oi6Y6gFLAYN2+dn3M265VPrjqEWsBgkFSyGDtL722d43SV/W3eMF7n0Lf/B9sUn1B1DTWYxdJLdvcz8ybN1p3iRM6Y8x9cu+jK733Jc3VHURBZDhxm3cxfrdrbXlY6vnziJZw7rhoi6o6hJLIYOk2vW8oF/Pq/uGIXbL/4a42cMeN9fdSCLQVLBYuhAh9z2NKevP73uGIWNl86oO4KaxGLoQLvv/zn8zWTO7mmvG6fcceLX6o6gJrEYOlTv+o2sfWpm3TE0RlkMHWzKwof4wKaT6c3ddUfRGGMxdLgnTvw1izb8Zd0xNMZYDGNA7zuf4bV3vL/uGBpDLIYxYPf27bxq2a848kd/VXcUjREWwxixa8vj/PG5WzjmzvfWHUVjgMUwhux6fCuHvu8RztxwWi0nJDf3+pXoscJiGGN2P/ccz//ZVpY88hejvu1PnXDmqG9TrWExjFFPnPQb3njXe/jwL99SdxR1IIthrMpk2jseZOtfz+TIb3+Ef9j2upZu7shbP8zu7c+0dBsaPRbDGNf7wIMc8cnb+enfzeeK37yiJduY84OzOOrcx9j9bHvdK0LDZzG8RIz7n//j+jPfzNqdzzd1vUff/n5e88lN9G7d1tT1ql4Ww0tI7/qNfGLeKSw84wOs3fk823qH/zf82p3Pc+q6dzJ78QZ6n/xVE1OqHYyvO4BG1+7nnoPV93Pu4W/i8Y+fyMIP3cZ50/+XQ7qmNLyOzz5xNHccN4lxux8lW5hV9RmyGCLiSuAdwLbMfG01Ng24Bjgc2AS8JzOfrn53PnAW0Av8fWb+oCXJNWKvuOQ27r4ETrj0PLoPbfzdw5xztsJuDx3Gssjcd+dHxFuBZ4B/61cMXwCeysyLImI5MDUzPx0R84CrgQXAocDNwFGZ2buvbRwU0/L4OGXkfxpJg7o5r1uTmQ39r8qHPMeQmT8GntpreBFwVfX8KuDMfuMrM3NHZj4MbKSvJCR1kOGefJyRmVsAqsc9dwGdBTzab15PNSapgzT75ONAX5Yf8FglIpYCSwEmMbnJMSSNxHDfMWyNiJkA1eOeM1E9wOx+8w4DBvxfImfmisycn5nzJ9A9zBiSWmG4xbAKWFI9XwLc0G98cUR0R8QcYC5w58giShptjXxceTVwMjA9InqAzwIXAddGxFnAL4F3A2Tm2oi4FngA2AUsG+oTCUntZ8hiyMzB7vwx4OeLmXkhcOFIQkmql1+JllSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVLAZJBYtBUsFikFSwGCQVhiyGiLgyIrZFxP39xj4XEY9FxN3Vz+n9fnd+RGyMiPURcVqrgktqnUbeMXwTWDjA+CWZeWz1cyNARMwDFgNHV8tcFhFdzQoraXQMWQyZ+WPgqQbXtwhYmZk7MvNhYCOwYAT5JNVgJOcYzomIe6tDjanV2Czg0X5zeqqxQkQsjYjVEbH6BXaMIIakZhtuMXwVOAI4FtgCXFyNxwBzc6AVZOaKzJyfmfMn0D3MGJJaYVjFkJlbM7M3M3cDl/OHw4UeYHa/qYcBm0cWUdJoG1YxRMTMfi/fBez5xGIVsDgiuiNiDjAXuHNkESWNtvFDTYiIq4GTgekR0QN8Fjg5Io6l7zBhE3A2QGaujYhrgQeAXcCyzOxtSXJJLROZA54CGFUHxbQ8Pk6pO4Y0pt2c163JzPmNzPWbj5IKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIKQxZDRMyOiFsjYl1ErI2Ij1bj0yLipojYUD1O7bfM+RGxMSLWR8RprfwDSGq+Rt4x7ALOy8zXACcAyyJiHrAcuCUz5wK3VK+pfrcYOBpYCFwWEV2tCC+pNYYshszckpl3Vc+3A+uAWcAi4Kpq2lXAmdXzRcDKzNyRmQ8DG4EFTc4tqYX26xxDRBwOHAf8FJiRmVugrzyAQ6pps4BH+y3WU41J6hANF0NEHABcD3wsM3+7r6kDjOUA61saEasjYvUL7Gg0hqRR0FAxRMQE+krh25n5nWp4a0TMrH4/E9hWjfcAs/stfhiwee91ZuaKzJyfmfMn0D3c/JJaoJFPJQK4AliXmV/q96tVwJLq+RLghn7jiyOiOyLmAHOBO5sXWVKrjW9gzknAB4H7IuLuauwC4CLg2og4C/gl8G6AzFwbEdcCD9D3icayzOxtdnBJrTNkMWTmTxj4vAHAKYMscyFw4QhySaqR33yUVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBUsBkkFi0FSwWKQVLAYJBWGLIaImB0Rt0bEuohYGxEfrcY/FxGPRcTd1c/p/ZY5PyI2RsT6iDitlX8ASc03voE5u4DzMvOuiDgQWBMRN1W/uyQzv9h/ckTMAxYDRwOHAjdHxFGZ2dvM4JJaZ8h3DJm5JTPvqp5vB9YBs/axyCJgZWbuyMyHgY3AgmaElTQ69uscQ0QcDhwH/LQaOici7o2IKyNiajU2C3i032I9DFAkEbE0IlZHxOoX2LH/ySW1TMPFEBEHANcDH8vM3wJfBY4AjgW2ABfvmTrA4lkMZK7IzPmZOX8C3fubW1ILNVQMETGBvlL4dmZ+ByAzt2Zmb2buBi7nD4cLPcDsfosfBmxuXmRJrdbIpxIBXAGsy8wv9Ruf2W/au4D7q+ergMUR0R0Rc4C5wJ3Niyyp1Rr5VOIk4IPAfRFxdzV2AfDeiDiWvsOETcDZAJm5NiKuBR6g7xONZX4iIXWWyCwO/0c/RMQTwLPAk3VnacB0OiMndE7WTskJnZN1oJyvysyDG1m4LYoBICJWZ+b8unMMpVNyQudk7ZSc0DlZR5rTr0RLKlgMkgrtVAwr6g7QoE7JCZ2TtVNyQudkHVHOtjnHIKl9tNM7BkltovZiiIiF1eXZGyNied159hYRmyLivurS8tXV2LSIuCkiNlSPU4daTwtyXRkR2yLi/n5jg+aq81L4QbK23WX7+7jFQFvt11G5FUJm1vYDdAG/AF4NTATuAebVmWmAjJuA6XuNfQFYXj1fDny+hlxvBd4A3D9ULmBetW+7gTnVPu+qOevngE8MMLe2rMBM4A3V8wOBB6s8bbVf95Gzafu07ncMC4CNmflQZu4EVtJ32Xa7WwRcVT2/CjhztANk5o+Bp/YaHixXrZfCD5J1MLVlzcFvMdBW+3UfOQez3znrLoaGLtGuWQI/jIg1EbG0GpuRmVug7x8ScEht6V5ssFztup+Hfdl+q+11i4G23a/NvBVCf3UXQ0OXaNfspMx8A/B2YFlEvLXuQMPQjvt5RJftt9IAtxgYdOoAY6OWtdm3Quiv7mJo+0u0M3Nz9bgN+C59b8G27rm6tHrcVl/CFxksV9vt52zTy/YHusUAbbhfW30rhLqL4WfA3IiYExET6btX5KqaM/1eREyp7nNJREwB3kbf5eWrgCXVtCXADfUkLAyWq+0uhW/Hy/YHu8UAbbZfR+VWCKNxtneIM6yn03dW9RfAZ+rOs1e2V9N3NvceYO2efMAfAbcAG6rHaTVku5q+t4sv0Pc3wln7ygV8ptrH64G3t0HWbwH3AfdW/+LOrDsr8Gb63mLfC9xd/Zzebvt1Hzmbtk/95qOkQt2HEpLakMUgqWAxSCpYDJIKFoOkgsUgqWAxSCpYDJIK/w/AEZImjwWu/QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "from scipy.interpolate import interpn\n", + "\n", + "img_shape = edt_1.shape\n", + "values = np.stack([edt_1, edt_2], axis=0)\n", + "points = tuple([np.arange(i) for i in values.shape])\n", + "xi = xi_coords(img_shape, percent=0.5)\n", + "interpolated = interpn(points, values, xi, method='linear')\n", + "result = np.reshape(interpolated, shape) > 0\n", + "\n", + "print(result.shape)\n", + "print(result)\n", + "plt.imshow(result)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 203, + "id": "3fda781b", + "metadata": {}, + "outputs": [], + "source": [ + "outputs = []\n", + "\n", + "percentages = np.arange(0, 1.1, 0.1)\n", + "for percent in percentages:\n", + " xi = xi_coords(img_shape, percent=percent)\n", + " interpolated = interpn(points, values, xi, method='linear')\n", + " out = np.reshape(interpolated, shape) > 0\n", + " outputs.append(out)\n", + "outputs = np.array(outputs)" + ] + }, + { + "cell_type": "markdown", + "id": "ed67a781", + "metadata": {}, + "source": [ + "## Visualization" + ] + }, + { + "cell_type": "code", + "execution_count": 205, + "id": "92e7849e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 205, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "viewer.add_labels(outputs)" + ] + }, + { + "cell_type": "code", + "execution_count": 206, + "id": "de1f65d7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 206, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "viewer.add_labels(image_1)\n", + "viewer.add_labels(image_2)" + ] + }, + { + "cell_type": "markdown", + "id": "b0afc2e9", + "metadata": {}, + "source": [ + "## Test cases\n", + "* a simple 3D test case (interpolating 2D slices)\n", + " * assert when percent = 0, the result exactly matches the input image_1\n", + " * assert when percent = 1.0, the result exactly matches the input image_2\n", + " * check results against a previous computation\n", + "* A test case with multiple label id values\n", + "* A 4D test case (interpolate across time points, from two labelled 3d blobs)\n", + "* Check code is robust as to which order image_2 and image_2 are given (eg: user can be scrolling forwards OR backwards through the slices)\n", + "* A test case when the structure is branching (one input blob turns into two output blobs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1474cb7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}