{ "cells": [ { "cell_type": "markdown", "id": "388df18a", "metadata": {}, "source": [ "# Methods and Parameters for Feature Detection: Part 2\n", "\n", "In this notebook, we will contninue to look in detail at tobac's feature detection and examine the remaining parameters.\n", "\n", "We will treat:\n", "\n", "- [Object Erosion Parameter](#Object-Erosion-Parameter-n_erosion_threshold)\n", "- [Minimum Object Pair Distance](#Minimum-Object-Pair-Distance-min_distance)\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "7451e92f", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:28.309652Z", "iopub.status.busy": "2026-02-02T20:12:28.309502Z", "iopub.status.idle": "2026-02-02T20:12:29.604538Z", "shell.execute_reply": "2026-02-02T20:12:29.604044Z" } }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", "\n", "%matplotlib inline\n", "\n", "import seaborn as sns\n", "\n", "sns.set_context(\"talk\")\n", "\n", "\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": 2, "id": "da405db0", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:29.606185Z", "iopub.status.busy": "2026-02-02T20:12:29.606035Z", "iopub.status.idle": "2026-02-02T20:12:30.364714Z", "shell.execute_reply": "2026-02-02T20:12:30.364091Z" } }, "outputs": [], "source": [ "import tobac\n", "import tobac.testing" ] }, { "cell_type": "markdown", "id": "58300a2b", "metadata": {}, "source": [ "## Object Erosion Parameter `n_erosion_threshold`" ] }, { "cell_type": "markdown", "id": "e25104dc", "metadata": {}, "source": [ "To understand this parameter we have to look at one varibale of the feature-Datasets we did not mention so far: *num*\n", "\n", "The value of *num* for a specific feature tells us the number of datapoints exceeding the threshold. *n_erosion_threshold* reduces this number by [eroding](https://en.wikipedia.org/wiki/Erosion_%28morphology%29) the mask of the feature on its boundary. Supose we are looking at the gaussian data again and we set a treshold of 0.5. The resulting mask of our feature will look like this:" ] }, { "cell_type": "code", "execution_count": 3, "id": "7cb13774", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.366475Z", "iopub.status.busy": "2026-02-02T20:12:30.366267Z", "iopub.status.idle": "2026-02-02T20:12:30.422271Z", "shell.execute_reply": "2026-02-02T20:12:30.421806Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAKnCAYAAABOG2+0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJ8ZJREFUeJzt3Q2QVeV9+PFneUeIKSggggqCoIEEAWuwSgjEiIwZUEln0mQ005CGNKnRlgwwmVqYTGqDpI0vmQbIMCYpNp1ImQCpDRVZR2mkBlCmYwhvxaCIaAQN7xC4/3nOzN0/K7vL8uLub/d+PjObc3PPubuHfXbly7nnPKeqVCqVEgAABNOmuXcAAADqIlQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACKlZQrW6ujp96lOfSj169EidO3dOV199dbr//vvTgQMHmmN3AAAIqKpUKpWa8gs++uij6d577035y/bt27eI1V//+tfpyJEj6ZprrkmrV69O3bt3P29fr1+/funNN99MnTp1Sv379z9vnxcAgLOzffv2dPjw4dSzZ8/0yiuv1L9hqQmtXbu21KZNm1JVVVVp/vz5pRMnThTP79y5szRy5MgczKU777zzvH7Nzp07F5/Xh++BnwE/A34G/Az4GfAz4Gcghfoe5E5rSJMeUb399tvT0qVL0913351+9KMf1Vq3ZcuW4hSAEydOpA0bNqSPfOQj5+Vr5qOze/fuTW1S29QlfeC8fE4AAM7egbQvnUjHU7du3dKePXvq3a5daiL79+9Pv/jFL4rHX/rSl05Zf9VVV6Vx48allStXpieeeOK8hWp+uz+Hao7Uj1bdfF4+JwAAZ+9/SivTvvTOaU/LbLKLqV588cXiPNSOHTum66+/vs5tRo8eXSzXrFnTVLsFAEBQTRaqmzdvLpaXX355at++fZ3bDBgwoFhu2rSpqXYLAICgmuyt//L5Bw1d0V9el9+qb8j8+fPTggULGvV1N27ceEb7CQBAhYVqnoIg69ChQ73b5NMCskOHDjX4uXbt2pXWr19/nvcQAICKDNU8j2l29OjRerfJ57Bm+SYADendu3caMWJEo4+oni58AQCo4FDN0w9kDU1BUF5X3rY+U6dOLT4aY+TIkY6+AgC0QE12MdWgQYOK5Y4dO9KxY8fq3Gbbtm21tgUAoHI1Wajmt+rz+an57f0XXnihzm2ee+65YnnDDTc01W4BAFDpodq1a9c0fvz44nFdV+znO1OtWrWqePzpT3+6qXYLAIBKD9Xs/vvvT1VVVelf/uVfilgt3701X8X/Z3/2Z8XtU/NtVocNG9aUuwUAQKWH6h//8R+nf/qnfyoe54uhrrjiiuKUgHz7rHXr1qXBgwenH/zgB025SwAABNWkoZrdd9996amnnkoTJkxIBw4cSL/+9a+LYP3GN76R1q5dmy6++OKm3iUAACp5eqqTfeITnyg+AAAgzBFVAABoDKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFAKD1hOobb7yRFi1alO6999504403pgsuuCBVVVWl66677rSvPXbsWJo7d24aNmxY6tKlS+revXsaN25cWrJkydnsCgAArVS7s3nRv/3bv6W//uu/PuPXHT58OH3yk59Mq1evTm3btk1DhgxJBw4cSNXV1cXHjBkz0re//e2z2SUAAFqZszqieuGFF6abb745zZw5My1evDg98MADjXpdDtEcqf37908vv/xy2rBhQ9q6dWtaunRp6tixY5ozZ05avnz52ewSAACtzFmF6he+8IX01FNPpX/4h39IkydPTr179z7ta3bv3p3mzZtXPF64cGEaPHhwzbqJEyem6dOnF49nz559NrsEAEAr02QXUy1btiwdPXo0DRw4MI0dO/aU9VOnTi2W69evT9u2bWuq3QIAoNJDdc2aNcVy9OjRda7v06dPcUrAydsCAFC5mixUN2/eXCzzEdX6DBgwoFhu2rSpqXYLAIDWdNX/2dizZ0+xzNNR1ae8bu/evQ1+rvnz56cFCxY06utu3LjxjPYTAIAKC9U8NVXWoUOHerfJV/5nhw4davBz7dq1qziXFQCA1qvJQrVTp07FMl9QVZ8jR44Uy86dOzf4ufIsAyNGjGj0EdXThS8AAPE0Wah269at1ikAdSmvK29bnzxDQHmWgNMZOXKko68AAC1Qk11MNWjQoGKZJ/ivT3laqvK2AABUriYL1VGjRhXLfGequuzcuTNt37691rYAAFSuJgvVSZMmpfbt26ctW7ak6urqOq/kz4YPH97gFFYAAFSGJgvVXr161ZxXOmXKlFpzpS5fvjw9+OCDxeNZs2Y11S4BANDaLqZ69dVXiyOf771af8OGDeniiy+ueX769OnFR1mO0XXr1qXnn38+DRkyJA0dOjTt37+/5tzUadOmFUdeAQDgrEL1+PHj6e233z7l+T/84Q+1nj948GCt9XnaqWeeeSY99NBDadGiRcXdqvK8qmPGjEn33HNPmjx5shEBAODsQ7Vfv36pVCqdzUuLMH3vkVYAAGi2c1QBAOBMCFUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACCkds29AwBNZcXrL/lmBzD+0mubexeAFsIRVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACKldc+8AQLbi9Zd8IypEU4z1+Euvfd+/BvD+c0QVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCatfcOwC0Tk1xP3c4Xz9/4y+91jcTAnJEFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACE1K65dwCIb8XrLzX3LkC4n/Hxl177vuwL8P85ogoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACG1a+4dAFrGfc2Bc/s9Gn/ptb6FcIYcUQUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQCgdYRqqVRKv/zlL9PMmTPTTTfdlC666KLUvn371KNHj3TLLbekxx9/vNimPseOHUtz585Nw4YNS126dEndu3dP48aNS0uWLDnXPwsAAJV8C9VVq1alm2++ueb/X3nllal///5p+/bt6amnnio+fvKTn6R///d/Tx07dqz12sOHD6dPfvKTafXq1alt27ZpyJAh6cCBA6m6urr4mDFjRvr2t799fv5kAABU3hHVHKYPP/xw2r17d9q2bVtau3Ztevvtt9OPf/zjIk7/4z/+I82aNeuU1+YQzZGaX//yyy+nDRs2pK1bt6alS5cWr5szZ05avnz5+fqzAQBQSaF6/fXXp02bNqWvfe1rqWfPnrXW3XXXXenv/u7visc/+MEP0okTJ2rW5aidN29e8XjhwoVp8ODBNesmTpyYpk+fXjyePXv22f9pAACo3FC98MILi3NS6zNhwoRiuWfPnvTWW2/VPL9s2bJ09OjRNHDgwDR27NhTXjd16tRiuX79+uIoLQAAle28X/Wfz0Mt69y5c83jNWvWFMvRo0fX+bo+ffoUpwScvC0AAJXrvIdqvpAqy1f156OvZZs3by6W+YhqfQYMGFAs86kFAABUtjO+6r8h+W378nmoefqqk+VTAbI8HVV9yuv27t3b4NeZP39+WrBgQaP2aePGjY3aDgCAVhqq+WKpO+64o5gnNS8/85nP1HlKQIcOHer9HOXprA4dOtTg19q1a1cRxQAAtF7nJVTffffd4iKqHTt2pJEjR6Yf/vCHp2zTqVOnYpkvqKrPkSNHTjm3tS69e/dOI0aMaPQR1dOFLwAArTBU9+/fn2699db04osvFhP4r1ixota5qWXdunWrdQpAXcrrytvWJ88QUJ4l4HRyODv6CgBQYRdTHTx4MN12223FVfqDBg1KK1euLG6pWpe8PssT/NenPC1VeVsAACrXWYdqPud00qRJ6dlnn039+vVLTz/9dLrkkkvq3X7UqFHFMt+Zqi47d+4sbsN68rYAAFSuswrVfMHU5MmTiyOoffv2TatWrSqWDclRm28UsGXLllRdXV3nlfzZ8OHDG5zCCgCAynDGoXr8+PH0uc99Lj355JPFEdQcqeWJ+hvSq1evmvNKp0yZUmuu1OXLl6cHH3yweDxr1qwz3SUAAFqhM76Y6qc//Wl64oknaq7k//M///N6t3300UeLI6RlOUbXrVuXnn/++eLCq6FDhxYXY5XPTZ02bVpx5BU4Mytef8m3DFrZ7+n4S6993/YFWm2olqeQyl555ZXio6Fpq06Wp5165pln0kMPPZQWLVpU3K0qz6s6ZsyYdM899xSnEwAAQFZVKpVKrflbUZ6e6gPpj9JHq25u7t2B94UjqtD6OKJKa/Y/pZVpX3qnmBc/v9v+vkxPBQAA7xehCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAIqV1z7wBw7sZfeu0Zbb/i9Zd82yH47yngiCoAAEF56x8AgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCA1hOqy5cvT1/96lfTqFGjUt++fVOnTp1S165d09ChQ9N9992Xfvvb39b72mPHjqW5c+emYcOGpS5duqTu3buncePGpSVLlpzLnwMAgFbmrEL1H//xH9M///M/p/Xr16e2bdumD3/4w6lHjx5p48aN6eGHH04f+tCH0n/913+d8rrDhw8XUTp9+vT08ssvp4EDBxahWl1dnSZPnpxmzpx5Pv5MAABUaqh+4QtfSCtXrkz79u0rjp7+6le/Stu3b0+bN29OH/vYx9LBgwfT5z73uXTgwIFar5sxY0ZavXp16t+/fxGqGzZsSFu3bk1Lly5NHTt2THPmzCmO1gIAwFmF6t13350+8YlPFHF5sgEDBqSf/vSnxePf/e536dlnn61Zt3v37jRv3rzi8cKFC9PgwYNr1k2cOLE4yprNnj3bqAAAcP4vpurVq1fxdn6Wj6yWLVu2LB09erR4u3/s2LGnvG7q1KnFMp9OsG3bNkMDAFDhznuo5vNU9+zZk9q0aZOGDx9e8/yaNWuK5ejRo+t8XZ8+fYpTAk7eFgCAytXufHySUqmU3nrrreL803weavb1r389XXnllTXb5PNXs3xEtT751IF8ruumTZsa/Hrz589PCxYsaHQ4AwBQYaG6aNGidNddd9V67uqrr06PP/54+uxnP1vr+XyUNSufFlCX8rq9e/c2+HV37dpVnCIAAEDrdU6h2rNnz3TjjTemEydOpNdeey3t3LmzOHKaQzVf/Z/nWD15aqqsQ4cO9X6+8sVZhw4davDr9u7dO40YMaLRR1RP9/kAAGhloXrLLbcUH2X/93//l6ZNm5Z+9rOfFTcDyFNQffCDHyzW5ZsCZPmCqvocOXKkWHbu3LnBr5svvCpffHU6I0eOdPQVAKDSL6bK56QuXrw4DRkypDi6+r3vfa9mXbdu3WqdAlCX8rrytgAAVK7zftV/vlPVhAkTisdr166teX7QoEHFMk/wX5/ytFTlbQEAqFznPVSzY8eOFct87mpZPhUgyzMD1CUfgc1X/J+8LQAAleu8h2o+B/XnP/958fjkeVQnTZqU2rdvn7Zs2ZKqq6vrnHKq/JqGprACAKAynPHFVPnt/KVLlxa3Ub3qqqtqrctX/N9zzz3FW/hdu3ZNf/EXf1HrjlX5Aqh83uqUKVPSf/7nf9bcRnX58uXpwQcfLB7PmjXr3P9UQIPGX3rtGX2HVrz+ku8onOPvEdAEobp///70rW99q/jo0aNHuuyyy4ojpXlu0x07dtTMh/rEE08Ud5s6WY7RdevWpeeff7644Gro0KHF5yufm5pnDMhHXgEA4IxDddiwYemRRx5JzzzzTPrf//3f4uKogwcPFtNQ3XTTTenWW28tjpxefPHFp7w2TzuVX/fQQw8VNwvIR2DzvKpjxowpjsROnjzZiAAAUKgq5fuftmLleVQ/kP4ofbTq5ubeHWiRvPUPp/LWP5y9/ymtTPvSO8UNnPK77U161T8AAJwroQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCANA67kwFVJ6zmdjcTQJoSUzeDzE5ogoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQmrX3DsAtE7jL732jLZf8fpL79u+UHnO9OcPiMkRVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACKldc+8AQFPdm33F6y/5ZlfIWAOtgyOqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQUrvm3gGApuIe8wAtiyOqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgDQekP1ySefTFVVVcVHv3796t3u2LFjae7cuWnYsGGpS5cuqXv37mncuHFpyZIl52M3AABoRc45VPft25e+/OUvn3a7w4cPF1E6ffr09PLLL6eBAwcWoVpdXZ0mT56cZs6cea67AgBAK3LOoTpjxoz06quvpttvv/20261evTr179+/CNUNGzakrVu3pqVLl6aOHTumOXPmpOXLl5/r7gAA0EqcU6jm8Jw3b16644470qRJk+rdbvfu3cV22cKFC9PgwYNr1k2cOLE4yprNnj37XHYHAIBW5KxDNb+V/8UvfjF17do1Pfroow1uu2zZsnT06NHi7f6xY8eesn7q1KnFcv369Wnbtm1nu0sAALQiZx2q3/zmN9OmTZvSAw88kPr06dPgtmvWrCmWo0ePrnN9fn0+JeDkbQEAqGxnFaovvfRScfX+9ddfn77yla+cdvvNmzcXy3xEtT4DBgwoljl+AQCg3Zl+C44fP56mTJlSPF6wYEFq0+b0rbtnz55ima/yr0953d69e0/7+ebPn1987cbYuHFjo7YDAKCFh+p3vvOd4lzSfAFUng+1seezZh06dKh3m3zlf3bo0KHTfr5du3YV+wAAQOt1RqG6ZcuW4sr8fD7prFmzGv26Tp06Fct8QVV9jhw5Uiw7d+582s/Xu3fvNGLEiEYfUW1M/AIA0IJDNU/sn4+Ofv/7308XXHBBo1/XrVu3WqcA1KW8rrxtQ/IsAeWZAk5n5MiRjr4CALT2UF23bl1xm9TPf/7zp6wrH7XMk/9fcsklxeN8a9Q/+ZM/SYMGDUr//d//XUzwX5/ytFR5WwAAOONzVEulUjGBf31OnDhRs778Vv+oUaPSY489VtwgoC47d+5M27dvr9kWAADOaHqqd955pwjVuj5yiGZXXHFFzXMf//jHi+fyXavat29fnONaXV1d51X82fDhwxucwgoAgMpxTrdQbaxevXrVnFOap7Y6ea7U5cuXpwcffLB4fCYXaAEA0Lqd8Vv/ZyvHaD7H9fnnn09DhgxJQ4cOTfv37685N3XatGnFkVcAAGiyI6rlaaeeeeaZNGfOnPShD32ouFvV7373uzRmzJi0ePHiYn5WAAAoqyrlk0lbsfL0VB9If5Q+WnVzc+8OAEDF+5/SyrQvvVPMi5/fcW/2I6oAAHAmhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCABCSUAUAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIQhUAgJCEKgAAIVWVSqVSasW6d++e9u7dm9qktqlL+kBz7w4AQMU7kPalE+l46tatW9qzZ0+93492rf07dfjw4WKZvxn70jvNvTsAALyn0yo2VHv27JnefPPN1KlTp9S/f/+a5zdu3JgOHTqUOnfunK655ppm3UfeX8a6chjrymGsK4vxbn22b99eRGrutIoO1VdeeaXO50eOHJnWr19fROq6deuafL9oOsa6chjrymGsK4vxrlwupgIAICShCgBASEIVAICQhCoAACEJVQAAQhKqAACEJFQBAAhJqAIAEJJQBQAgJKEKAEBIrf4WqvX50pe+lHbt2pV69+7d3LvC+8xYVw5jXTmMdWUx3pWrqlQqlZp7JwAA4L289Q8AQEhCFQCAkIQqAAAhVVyoVldXp0996lOpR48eqXPnzunqq69O999/fzpw4EBz7xpn4I033kiLFi1K9957b7rxxhvTBRdckKqqqtJ111132tceO3YszZ07Nw0bNix16dIlde/ePY0bNy4tWbLEGASTT6H/5S9/mWbOnJluuummdNFFF6X27dsXv7+33HJLevzxx4tt6mOsW5bly5enr371q2nUqFGpb9++qVOnTqlr165p6NCh6b777ku//e1v632tsW75nnzyyeK/4/mjX79+9W5nrCtMqYI88sgjpaqqqvy3Wqlv376l4cOHlzp27Fj8/2uuuab09ttvN/cu0kjf/e53i3F778fIkSMbfN2hQ4dKN910U7Ft27ZtSx/5yEdKAwYMqHn9jBkzjEEgK1eurDW+V155ZTHG3bt3r3nutttuKx0+fPiU1xrrlmfMmDHFmLZv3750+eWXl6677rpSv379Sm3atCmev+CCC0orVqw45XXGuuX7/e9/X7rssstqfq+vuOKKOrcz1pWnYkJ17dq1xX/scqjOnz+/dOLEieL5nTt3Fn/x5V+MO++8s7l3k0ZauHBh6eabby7NnDmztHjx4tIDDzzQqFD92te+VmzXv3//0m9+85ua55cuXVrzj5Zly5YZhyCeeuqpYqwefvjh0u7du2ut+/GPf1wzZnX9A8NYtzw/+tGPin+cvPcfHlu3bi197GMfK8b64osvLu3fv7/WemPd8v3lX/5lMb633357g6FqrCtPxYTqpEmTih/+u++++5R1mzdvrvkX+4YNG5pl/zg3jz322GlD9Y033ih16NCh2G7VqlWnrL///vuLdSNGjDAcQbz77rulo0eP1rv+7//+74sxy0dYjx8/XvO8sW598piWj7Y9+eSTtZ73e92yPffcc8VBpDvuuKPmv+V1haqxrkwVcY7q/v370y9+8YuaSYPf66qrrirOUcyeeOKJJt8/msayZcvS0aNH08CBA9PYsWNPWT916tRiuX79+rRt2zbDEsCFF15YnJNanwkTJhTLPXv2pLfeeqvmeWPd+vTq1as4nzw7ePBgzfPGumU7fPhw+uIXv1ici/zoo482uK2xrkwVEaovvvhiOnLkSOrYsWO6/vrr69xm9OjRxXLNmjVNvHc0lfLYlsf6vfr06ZP69+9fa1vi/yVXli+OLDPWrc/GjRuLf5C0adMmDR8+vOZ5Y92yffOb30ybNm1KDzzwQPHf4IYY68pUEaG6efPmYnn55ZfXe3RmwIABxTL/wtC6fw7yEdX6+DloWX7yk58UyzyDQz76WmasW4d8etqbb75ZzMgxceLE4rmvf/3r6corr6zZxli3XC+99FIxA0s+gPSVr3zltNsb68rULlWA/K/wrPy2UV3K6/bu3dtk+0XT8nPQuuRTNObNm1c8ztNXncxYt2x56rm77rqr1nN5KsE8HdlnP/vZWs8b65bp+PHjacqUKcXjBQsWFEfKT8dYV6Y2lfT2YIcOHerdJp8WkB06dKjJ9oum5eeg9di9e3e64447ivkU8/Izn/lMrfXGumXr2bNnMT/yDTfckC677LIiYvLRtByqr732Wq1tjXXL9J3vfKf4x+bf/M3fFO+INIaxrkwVEap50ugsX0hTn3wO63vPc6N18XPQOrz77rvFRVQ7duxII0eOTD/84Q9P2cZYt2z5Zg6rV68ubvaQx3nLli3FW/95Qvh8M4D8M1BmrFuePJ6zZ88urgmYNWtWo19nrCtTRYRqt27dar1tUJfyuvK2tD5+DlrHDB633nprcYHkkCFD0ooVK2qdm1pmrFuXfE7q4sWLizHfuXNn+t73vlezzli3PF/+8peLo6Pf//73i7sKNpaxrkwVEaqDBg0qlvlf5vmtwrqUpyMqb0vrUx7brVu31ruNn4O48pREt912W3Hlbx7LlStXFrdUrYuxbn3atm1bMx3Z2rVra5431i3PunXritukfv7zn0+XXHJJrY98W+zs1VdfrXkuH1nPjHVlqohQHTFiRHF+an57/4UXXqhzm+eee65Y5nOiaJ3yW4ZZfkuxLvlIzfbt22ttSwz56MukSZPSs88+W9wD/Omnny7+AquPsW6dygcaTpw4UfOcsW65Mzrkc83f+/H73/++ZozLz5VP2zPWlakiQjVPJDx+/PiaqwvrOl9m1apVxeNPf/rTTb5/NI0cOnl6sjze1dXVp6yfP39+scxzNDY0hRVNHyeTJ08ujqD27du3+F3Ny4YY69Ynx8rPf/7z4vHJ86ga65bnnXfeKUK1ro/HHnus2OaKK66oee7jH/948ZyxrlClCvHCCy8Ut2jLH/Pnzy+dOHGieP71118vbrtZvscwrfcWqtlf/dVfFdvl+8f/5je/qXl+2bJlNfeN/9nPftYEe0xj/OEPfyj96Z/+aTEul1xySXG748Yy1i3Lr371q9Lf/u3f1jnGmzZtKt1yyy3Fz0HXrl1Lr732Wq31xrr1aOgWqpmxrjwVE6rZd7/73SJU8y/BZZddVho+fHhNnAwePLj01ltvNfcu0kg7duwoXXTRRTUf+S+vPI7t2rWr9fycOXNqve7gwYOlG264odi2bdu2pWHDhpUGDBhQcw/xadOmGYNA/vVf/7VmbPr161e68cYb6/1Yv359rdca65alurq6Zqx79OhRGjFiROmjH/1o6fLLL695vnv37qWnn376lNca68oJVWNdeSoqVLOVK1eWJkyYUPwHL0fqoEGDSt/4xjdK+/bta+5d4wxs37695i+vhj5mzZp1ymuPHDlSBOyHP/zhUufOnUsf/OAHS2PGjCktXrzYGAT9S6sxHzl03stYtxx79uwpPfLII6U777yzdNVVV5UuvPDCmn943nTTTaVvfetbDR5MMNaVEaqZsa4sVfl/mvv0AwAAqMiLqQAAaHmEKgAAIQlVAABCEqoAAIQkVAEACEmoAgAQklAFACAkoQoAQEhCFQCAkIQqAAAhCVUAAEISqgAAhCRUAQAISagCAJAi+n+SeAQhS4icywAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = np.linspace(-2, 2)\n", "y = np.linspace(-2, 2)\n", "xx, yy = np.meshgrid(x, y)\n", "\n", "exp = np.exp(-(xx**2 + yy**2))\n", "\n", "gaussian_data = np.expand_dims(exp, axis=0)\n", "threshold = 0.5\n", "\n", "mask = gaussian_data > threshold\n", "mask = mask[0]\n", "\n", "plt.figure(figsize=(8, 8))\n", "plt.imshow(mask)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d7e0e85a", "metadata": {}, "source": [ "The erosion algorithm used by tobac is imported from *skimage.morphology*:" ] }, { "cell_type": "code", "execution_count": 4, "id": "ae4d979f", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.423657Z", "iopub.status.busy": "2026-02-02T20:12:30.423557Z", "iopub.status.idle": "2026-02-02T20:12:30.501988Z", "shell.execute_reply": "2026-02-02T20:12:30.501367Z" } }, "outputs": [], "source": [ "from skimage.morphology import binary_erosion" ] }, { "cell_type": "markdown", "id": "09500e57", "metadata": {}, "source": [ "Applying this algorithm requires a quadratic matrix. The size of this matrix is provided by the *n_erosion_threshold* parameter. For a quick demonstration we can create the matrix by hand and apply the erosion for different values:" ] }, { "cell_type": "code", "execution_count": 5, "id": "43921591", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.503690Z", "iopub.status.busy": "2026-02-02T20:12:30.503465Z", "iopub.status.idle": "2026-02-02T20:12:30.749224Z", "shell.execute_reply": "2026-02-02T20:12:30.748727Z" } }, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, '$\\\\mathtt{n\\\\_erosion\\\\_threshold} = 10$')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABH8AAAGUCAYAAACskI5kAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMqpJREFUeJzt3QuclXWdP/DfgNzES6CIBAoIgYpGMKZ4QRSvaAlGbZYv7aKl22bauqtmEXb3lnlpS1Bfmou1m8aGtO6aBKxRoskkr9aIW5hGSCao3EF4/q/nt/8zLwZmhplh5szM73m/X6/jczzP85zzO+c5fM453/k9v19FlmVZAAAAACBJHVq7AQAAAAC0HMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8Yc247TTTgsVFRUhNS+99FJ8Xh//+MdDe9PSbf/ud78b73/nS/6YjTVgwIB4Adof2d/2yH6gOcj3tke+F5viD+wlhYemO/7448PkyZPjZfjw4e3yvfib3/wmnHfeeaFHjx6he/fu8Tn98Ic/bO1mAS1M9hc3+/Njv+sfLkqXK6+8srWbB+wl+V7cfJ82bVq44oorwnHHHRe6dOkSc/2hhx5K5rfAPq3dACh5+OGHw8aNG70gBZKHY34p/SVi4cKFoT2ZO3duOOecc0Lnzp3DRRddFA488MAwffr0cPHFF8fnc+ONN7Z2E6HNk/3F096zP5fn/TXXXLPb7fkPBuD/yPfiae/5/qUvfSn86U9/CgcffHDo06dPvJ7SbwHFH9qMww8/vLWbAA329ttvh8svvzz+ReDpp58OI0aMiLfnf+k48cQT4/JDH/pQeNe73uVVhXrIftqjd7zjHeGmm25q7WZAmybfaW/uv//++N29f//+4eabbw5f+MIXkvot4LSvMp9bOW/evHDKKaeEfffdNxx66KGxGrhjx469foyf/exnYezYsbHamN/3qFGjwowZM+rcPm9LaXyV/E1+zDHHhG7duoXevXvHrm4ledvuuOOOuL5r166xO9v5558fu7fV5dFHH43PsVevXrHr29ChQ+M/jGXLlu22bf7Fadcu03VpTFta+jUvncOcX/KKcH7Z+Tnk62vTmLY09Bg15T3QmGPU2LY35T3TUOvWrQtXX311rMTn7Tj11FNDVVVVaA2zZ88Oy5cvDx/96Eerwz63//77h0mTJsUPhAcffLBV2kbbIPtlf4nsTyf7ISff5XuJfE8r388888xY+En2t0BGi1uxYkWWv9QnnHBCdsABB2QXX3xx9rnPfS575zvfGW+/66679ur+b7755ng/hxxySHb55ZdnV111VTZ48OB424MPPljrPh/72Mfi+s9//vPZvvvum11yySXZ9ddfn33gAx/IjjzyyOrtLrvssrhdfts///M/Z5/85Cezbt26ZV27ds3mzZu32/1+73vfi9sPGDAg++xnPxv3ufDCC7P99tsv+9GPfrTb9nPmzMkmT54cL/3794/71qUxbWnp1zx/XUvtPvDAA+Ol9P/5ZefXvaltaegxaux7oDHHqCltb+x7Ztfnmz9mbbZv356NGTMmbnPqqadmN9xwQ3bWWWdlPXv2zHr06BHfP+X0hS98Ibaltvf1mjVr4rqTTjqprG2ibZH9sl/2p5f9ufwxDz300Oyhhx7KvvGNb8TP1RdeeKHs7aD1yHf5Lt/TzPedfetb36r393R7/C2g+FPGD4j88sQTT1TfvnTp0qxDhw7Zcccd1+T7rqqqiveR/9B+/fXXq2/fuHFjNmLEiFiUWLduXZ3/GA866KBs8eLFNda9+uqrcfncc8/Fbd797ndnmzZtql4/e/bseHt+/7uqrKyMQbB27doat7/11lvZqlWr6n0upX/8tWlsW1ryNd9VHkz1hVNT29KQY9SU90BjjlFj296U90xDPyD+/d//Pa6fMGFCrcWmhnxA7Fyga8hl19doZx/84Afj4z7//PO1rj/44IOzXr167bFNpEv2y37Zn17250p/rNr1cu6552avvfbaHttD+yff5bt8TzPfG1P8aY+/BRR/yvgBceyxx+627uijj449LprqyiuvjPf9H//xH7utu+++++pcV/rHOGnSpDrvO6++5tvcf//9u61773vfG9ctX768xu3vec97su7du2fr169v9HOpr/jT2La05Gve1OJPY9vSkGPUlPdAY45RY9velPdMQz8gPvzhD8f1eW+xnS1atKjBHxC1fVmv71JXW3L5XybybfJCWG2OOOKIrHPnzntsE+mS/Q0j+2uS/W07+3Nf+cpXsrlz58ZCT/6Hk/nz52fjxo2L+5544onZjh07Gvjup72S7w0j32uS720/3xtT/GmPvwUM+FxG+bgqu8pHEv/973/f5PssjaOSDzL1wgsv1FhXGr+lvnFcTj/99DrXLVq0KC5rm6YvP68xf+x8myOOOKL69okTJ8ZzHN/73veGj3zkI3Hcmfx6PjDi3mhKW1rqNW+qpralvmPUlPdAU45RQ9ve1OPUEKX7PvbYY2vcfuSRR8apGBvi/z4joLxkf9PJ/rrJ/tbL/i9/+cs1/v+EE06I4+6NGTMmjv/xxBNPxLHuSJ98bzr5Xjf57rt9S1H8KaN8YN1d1TfAcUO88cYbcfmd73ynzm02bNhQ57p+/frVuW79+vVx2bNnz93WHXTQQdWDdO0sHwQ4LyI88MADcYTz/AtXp06d4kjn9957bxwAqyma0paWes2bqqltqe8YNeU90JRj1NC2N/U4NUTpOeQDSO+qtsdrafnA2rk333yz1vVvvfVW9TYUm+yX/buS/e03++vSoUOH8IlPfCIWf371q18p/hSEfJfvu5Lv6eV7Sr8FFH/auVJvjfzN1ZTCSv6jvy777bdfXK5Zs2a3nhqvv/56XO76mPmXn89+9rPxkv9D+J//+Z9w++23hx/+8IfxH/E999zT6DY2tS2pqO8YNeU90FLHqKWPU+m+165dW11IKskfL5+BbE8aOy3vNddcU2ePqNK0jUuXLg2VlZU11uVt/Nvf/hZOOumkRj0eNJTsl/2yv3Wyvz55j9jcxo0bG70vlMh3+S7f216+p/JbQPGnnTvuuOPCggULwvPPP7/H04MaKz+dJp8qfOHChfFxdlY6vSjfpi55pfOCCy4I5557bvxCNGfOnFZrS0vq2LFj2LJlS2iP74HmPEZ7e5xKxZ26vjQfffTR8X5/97vfhdNOO61Gt+GGvv5f+cpXQmN8/OMfr/MDIu/e/61vfSv8/Oc/DxdddFGNdfltpW2gJch+2S/7Wyf76/Pss8/G5YABAxq9L5TId/ku39tevqfyW6BDazeAvXP55ZfH7oXXXnttrC7uKu/VUd9pX/W58MIL4/Luu++ucapO/mZ+7rnn4rgugwYNqrHPk08+GbZv317jtpUrV8Yf9LV16WvJtpRLXjRZvXp19elXbf090FLHaG+P0+DBg+Ny/vz5ta7PxyrK3XHHHWHHjh3Vt3/7299ucPv+/yD3Db7U9wX+jDPOiL2b8h5TO4+1lD/vr33ta2GfffaJHzDQEmS/7Jf9rZP9+Vh3tX3e56d75Z9P+ThEH/jAB5r8bxvku3yX762T743VHn8L6PnTzuWV4a9//evhi1/8Yhx07rzzzgvvfOc7w8svvxx/bP/xj38Mq1atqvWc5D3JBzDM37APPfRQGDlyZBg3blwsLkyfPj1+ucl/4O/qwx/+cHys0aNHh/79+8cub/n2+Y/1vDixs5deeine987/v2v3vQkTJoT3vOc9TWpLueSvef5a54M75l/48uefH4O8R01bfA805hg11t4cp4svvjgOonn11VfHAM1PQcsr83n3zFz+2uY9m2bOnBm7UJ566qmxJ1De62lvi1ZNkQf6/fffH84555z4WuaDZx9wwAHxua5YsSIekyFDhpS9XRSD7Jf9sr91sv/HP/5xuPXWW+OX/vxHRP7Z9r//+7/xjxz5adX52HmHH3542dtFOuS7fJfvrZPvufy7fV7Mz+VnG5Rumzt3bvVv0/zSbn8LtPZ0Y0WaDjKf3q8xUyA2xn/9139l5557btazZ884pVw+Nd748eOzH/zgB9m2bdsaPa12ydtvv53ddtttcVrv/H4PPPDA+Dj5tKa1+Zd/+ZfsvPPOy/r165d16dIl69u3b3b++ednTz/99G7b5tP67Wk6vp2n1mtMW8rxmpds3rw5u/rqq+Nz7dChQ7zv/DH2ti0NPUaNfQ805hg1pe2Nfc/sLN/mpJNOyrp161brFI/5lLpXXXVV1rt376xr167ZySefnC1YsCBu15DpIFvCs88+G59f/jzzdh933HHZtGnTWqUttC2yX/bL/vSyP5/i/e/+7u+ywYMHZ/vvv3/WqVOn+Hl60UUXxc8DikG+y3f5nl6+7/z7q67L5MmTs/b8W6Ai/09rF6AAAAAAaBnG/AEAAABImOIPAAAAQMJapfiTTyf9vve9L/Tq1St069YtTv08adKkJs9KBUDbJ/sBikf2A7QNZR/z55577omz+eQP269fv1gAyqfN3LJlSzjqqKPi6Nr5LD9F1Jip5/JZJh544IEWbU9ReN2h5cl+GdTWyH5oeUXPfjnjdYfCFn/yKZmPP/74+AGQT4X5qU99KlRUVIS//OUvcVrsfH0+nfNPfvKTcjUJgBYm+wGKR/YDFLj4M2HChDBjxoxw6aWXhh/84Ac11i1dujSe/rVjx46wcOHC8O53v7tczQKgBcl+gOKR/QBtyz7leqD169eH//7v/47XP/3pT++2/l3velcYO3ZsmDVrVnj00UebrfiTd7f861//Grp27RoGDhzYLPcJ0NatWLEibN68ORxyyCHhpZdearV2yH6A8ily9vvODxTVigZmf9mKP7/97W/j+b1dunSJp37VZvTo0fFDYP78+c32uHnhZ9OmTfGydu3aZrtfgPYgz8DWJPsByq+I2e87P1B0f91D9pet+LNkyZK4PPzww0OnTp1q3WbQoEFxuXjx4mZ73LzHT1746RA6hu5h/2a7X4C2bENYF3aE7TEDW5PsByifIme/7/xAUW1oYPaXrfizZs2auKxvRP/Suj310JkyZUqYOnVqgx437/6Uyws/J1Sc2YgWA7Rfz2azwrrwRquf7ir7Acontez3nR+g+bK/bMWfUhGmc+fOdW6Tdw3N5T116rNq1apQVVXVzC0EoLnJfoDiaa7s950foPmUrfhT6oK0devWOrfJzw3OdevWrd776tOnTxg5cmSDHnfRokV7LCYB0DJkP0DxNFf2+84P0A6LPz169KjRDbQ2pXWlbetyxRVXxEtDVFZW6iUE0EpkP0DxNFf2+84P0Hw6hDIZMmRIXL788sth27ZttW6zfPnyGtsC0L7JfoDikf0ABS7+5Kdp5ef95l08n3vuuVq3+eUvfxmXJ554YrmaBUALkv0AxSP7AQpc/Nlvv/3COeecE6/XNlPX0qVLw+zZs+P1D37wg+VqFgAtSPYDFI/sByhw8Sc3adKkUFFREf71X/81FoCyLKseyf8jH/lI2LFjR5gwYUIYPnx4OZsFQAuS/QDFI/sBClz8ee973xvuuOOO6gHc+vfvH7uF5vPRL1iwIAwdOjTcd9995WwSAC1M9gMUj+wHKHDxJ3fNNdeEp556KowbNy5s2LAh/P73v49FoBtvvDE8//zz4eCDDy53kwBoYbIfoHhkP0ABp3rf2RlnnBEvABSH7AcoHtkPUNCePwAAAACUj+IPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACSsScWfV199NUybNi1cffXV4eSTTw777rtvqKioCMcdd9we9922bVu47bbbwvDhw0P37t1Dz549w9ixY8P06dOb0hQAykT2AxSL3AdIxz5N2enf/u3fwuc///lG77d58+Zw1llnhXnz5oWOHTuGYcOGhQ0bNoQ5c+bEy/XXXx9uvvnmpjQJgBYm+wGKRe4DFLznzwEHHBDOPPPMcMMNN4THHnssfPOb32zQfnlxJy/8DBw4MLz44oth4cKFYdmyZWHGjBmhS5cu4ZZbbgkzZ85sSpMAaGGyH6BY5D5AwYs/n/zkJ8NTTz0VvvWtb4WJEyeGPn367HGf1atXh3vvvTdef+CBB8LQoUOr111wwQXhuuuui9dvuummpjQJgBYm+wGKRe4DpKNsAz4//vjjYevWrWHw4MHh9NNP3239FVdcEZdVVVVh+fLl5WoWAC1I9gMUi9wHKHjxZ/78+XE5evToWtf37ds3ng6287YAtG+yH6BY5D5AwYs/S5Ysicu8509dBg0aFJeLFy8uV7MAaEGyH6BY5D5AQrN9NcWaNWviMp/avS6ldWvXrq33vqZMmRKmTp3aoMddtGhRo9oJQPOR/QDFIvcBCl78yad5z3Xu3LnObfIZv3KbNm2q975WrVoVxwYCoG2T/QDFIvcBCl786dq1a1zmgz7XZcuWLXHZrVu3eu8rn11s5MiRDe75s6diEgAtQ/YDFIvcByh48adHjx41uoLWprSutG1d8pnBSrOD7UllZaVeQgCtRPYDFIvcByj4gM9DhgyJy2XLltW5TWmK99K2ALRvsh+gWOQ+QMGLP6NGjYrLefPm1bp+5cqVYcWKFTW2BaB9k/0AxSL3AQpe/Bk/fnzo1KlTWLp0aZgzZ06tM3jlRowYUe908AC0H7IfoFjkPkDBiz+9e/euHqfnsssuC4sXL65eN3PmzHDrrbfG65MnTy5XkwBoYbIfoFjkPkBCAz6/8sorsYfOrrN0LVy4MBx88MHVt1933XXxUpIXeBYsWBCeeeaZMGzYsHDMMceE9evXV4/1c+2118a/FgDQ9sh+gGKR+wAFL/5s3749vP7667vd/vbbb9e4fePGjTXW51O4z507N9x5551h2rRpYcmSJaFz585hzJgx4aqrrgoTJ05sSnMAKAPZD1Asch+g4MWfAQMGhCzLmvSAebFn1x5BALR9sh+gWOQ+QDrKNuYPAAAAAOWn+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkLB9WrsB0Fqe/MsLSb3457zzPa3dBAAAANogPX8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYWb7Ihmpzd7V0s/f7GAAAADFoOcPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhpnqnXSn6dO7lei1NAw/QsLxsTrIXAGgpev4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwsz2RZtkVq+2+fqbiQZo79ry50tT2iaXAYCG0PMHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwU73TatrydLs0/piZbhhoS4ryGVPX85TJAMDO9PwBAAAASFijiz9ZloVf//rX4YYbbginnHJKOOigg0KnTp1Cr169wtlnnx0eeeSRuE1dtm3bFm677bYwfPjw0L1799CzZ88wduzYMH369L19LgC0ENkPUDyyH6DAp33Nnj07nHnmmdX/f8QRR4SBAweGFStWhKeeeipefvSjH4Wf/OQnoUuXLjX23bx5czjrrLPCvHnzQseOHcOwYcPChg0bwpw5c+Ll+uuvDzfffHPzPDMAmo3sByge2Q9Q8J4/ebHnrrvuCqtXrw7Lly8Pzz//fHj99dfDww8/HAs+//mf/xkmT5682755cScv/OT7v/jii2HhwoVh2bJlYcaMGXG/W265JcycObO5nhsAzUT2AxSP7AcocPHn+OOPD4sXLw6f+9znwiGHHFJj3SWXXBK+/OUvx+v33Xdf2LFjR/W6vFB07733xusPPPBAGDp0aPW6Cy64IFx33XXx+k033dT0ZwNAi5D9AMUj+wEKXPw54IAD4hg/dRk3blxcrlmzJrz22mvVtz/++ONh69atYfDgweH000/fbb8rrrgiLquqqmJvIgDaDtlPW57tqrZL0XldaA6yHyAdzT7bVz6uT0m3bt2qr8+fPz8uR48eXet+ffv2jaeD7bwtAO2D7AcoHtkPkPCAz3uSD/acy2fzyv9aULJkyZK4zHv+1GXQoEFx4Oj8tLL6TJkyJUydOrVB7Vm0aFEDWw5AU8l+gOJp6ez3nR+gjRZ/8lO2SuP65FPB7yw/DSyXT+1el9K6tWvX1vs4q1atio8FQOuT/QDFU47s950foA0Wf/IBnS+88MKwbdu2uLzoootq7RbauXPnOu+jNDX8pk2b6n2sPn36hJEjRza458+e7g+AppH9AMVTruz3nR+gjRV/3nzzzTjQ88svvxwqKyvDQw89tNs2Xbt2jct80Oe6bNmyZbexgmqTDw5dGiB6T/L26CUE0PxkP0DxlDP7fecHaEPFn/Xr14dzzz03/Pa3vw3Dhg0LTz75ZI1zfkt69OhRoxtobUrrStuSBrOuFPs4n/PO95S9LbQ82U85+Rwpz2spr9kT2Q/N+5kkd2k3s31t3LgxnH/++XF2riFDhoRZs2aFgw46qNZt8/W5ZcuW1Xl/pSneS9sC0PbIfoDikf0ABS3+5Ofyjh8/Pjz99NNhwIAB4Re/+EU49NBD69x+1KhRcTlv3rxa169cuTKO+L/ztgC0LbIfoHhkP0BBiz/54G4TJ06MPX369esXZs+eHZf1yQtFnTp1CkuXLg1z5sypdSrH3IgRI+qdFhKA1iH7AYpH9gMUtPizffv2cPHFF4cnnngi9vTJCz8DBw7c4369e/euHqT5sssuC4sXL65eN3PmzHDrrbfG65MnT25skwBoYbIfoHhkP0CBB3z+8Y9/HB599NHqkfw/8YlP1LntPffcE3vylOQFngULFoRnnnkmDg59zDHHxIHjSmP9XHvttbGHEABti+wHKB7ZD1Dg4k9pWsbcSy+9FC/1TQW5s3wqx7lz54Y777wzTJs2LSxZsiR07tw5jBkzJlx11VXxVDIA2h7ZD1A8sh8gHRVZlmUhYZWVlaGqqirsH94RTqg4s7WbU0im6C02U1i2jmezWWFdeCOMHDky9rgsGtmfFp8j5SGv278iZ7/cJ6XPHXlMS2T/Xk31DgAAAEDbpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhjZ7qHQAgxdlViq6u19+sM0Cq2urnTn3tksk0lZ4/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEmeodAACAJLXV6dyb+/mYAp490fMHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhij8AAAAACduntRsAAJA7553vqfWFePIvL3iBWvH1B0g129rj54uspqn0/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEme0LAACAwmmrs0ya0YuWoOcPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhpnqnsFMo0rxMSQm0Rr74LGm+1xKA8n7uyGTKSc8fAAAAgIQ1qfgzc+bM8A//8A9h1KhRoV+/fqFr165hv/32C8ccc0y45pprwp/+9Kc69922bVu47bbbwvDhw0P37t1Dz549w9ixY8P06dP35nkA0MJkP0DxyH6AAhd/vv3tb4fvfe97oaqqKnTs2DEce+yxoVevXmHRokXhrrvuCkcffXT4+c9/vtt+mzdvjoWe6667Lrz44oth8ODBsfgzZ86cMHHixHDDDTc0x3MCoAXIfoDikf0ABS7+fPKTnwyzZs0K69ati718fvOb34QVK1aEJUuWhFNPPTVs3LgxXHzxxWHDhg019rv++uvDvHnzwsCBA2PxZ+HChWHZsmVhxowZoUuXLuGWW26Jf10AoO2R/QDFI/sBClz8ufTSS8MZZ5wRCzY7GzRoUPjxj38cr//tb38LTz/9dPW61atXh3vvvTdef+CBB8LQoUOr111wwQWxN1DupptuatozAaBFyX6A4pH9AGlo9tm+evfuHU/lWrNmTewBVPL444+HrVu3xlO9Tj/99N32u+KKK8LXvva1eCrZ8uXLYyEJgPZB9tNazCjZuNcFmpPsp4jkK+1Vs8/2lY/7kxd+OnToEEaMGFF9+/z58+Ny9OjRte7Xt2/feDrYztsC0D7IfoDikf0ABev5k2VZeO211+J4Pvm4Prl/+qd/CkcccUT1Nvl4QLm8509d8t4++dhBixcvrvfxpkyZEqZOndrgDyUAmp/sByiecma/7/wAbaT4M23atHDJJZfUuO3II48MjzzySPjoRz9a4/a8N1AuPyWsLqV1a9eurfdxV61aFU8PA6D8ZD9A8bRG9vvOD9BGij+HHHJIOPnkk8OOHTvCn//857By5cpY6c8/BPJZv/r161djmvdc586d67y/0gDSmzZtqvdx+/TpE0aOHNngnj97uj8AGk72AxRPa2S/7/wAbaT4c/bZZ8dLyR//+Mdw7bXXhp/+9Kdh1KhRcTr3Aw88MK7r2rVrXOaDPtdly5YtcdmtW7d6HzcfHDq/NERlZaVeQgDNSPYDFE9rZL/v/ABtdLav/Fzfxx57LAwfPjx+AHz3u98NX/ziF+O6Hj161OgGWpvSutK2FHek/Cf/8kJZ20LDmN2A2sh+2qKizAIml2ktsh+g4LN9dezYMYwbNy5ef/7556tvHzJkSFwuW7aszn3zKd533haA9kH2AxSP7AcocPEnt23btrjMzwkuybuD5vKZAWqTnzecj/i/87YAtB+yH6B4ZD9AQYs/+bm9P/vZz+L1ESNGVN8+fvz40KlTp7B06dIwZ86cWqdyLO1T37SQALQ9sh+geGQ/QMLFn/xUrkmTJsUizq7yEf/f//73x9O39ttvv/CpT32qel3v3r2rB2m+7LLLwuLFi6vXzZw5M9x6663x+uTJk5v6XABoIbIfoHhkP0CBB3xev359+PrXvx4vvXr1Cocddljs0bNq1arw8ssvx2169uwZHn300dC3b98a++YFngULFoRnnnkmDBs2LBxzzDHx/kpj/eQzBuQ9hABoW2Q/QPHIfoACF3/ymbzuvvvuMHfu3PC73/0uDuC8cePGOLXjKaecEs4999zYw+fggw/ebd98Ksd8vzvvvDNMmzYt9hTq3LlzGDNmTLjqqqvCxIkTm+t5AdCMZD9A8ch+gHRUZFmWhYRVVlaGqqqqsH94Rzih4szWbg57KbUpetsqUwe3f89ms8K68EYYOXJk7HFZNLKf9vjZI3vZW0XOfrkPFNWzDcz+FpntCwAAAIC2QfEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCGj3VO7Sm+mZCMRNY872WAMhLACAdev4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABJmqncKO3V5alPDm7odAACA2uj5AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAkz2xeFZXYsAAAAikDPHwAAAICEKf4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhzVL8eeKJJ0JFRUW8DBgwoM7ttm3bFm677bYwfPjw0L1799CzZ88wduzYMH369OZoBgBlJPsBikf2AxS0+LNu3bpw5ZVX7nG7zZs3x0LPddddF1588cUwePDgWPyZM2dOmDhxYrjhhhv2tikAlInsByge2Q9Q4OLP9ddfH1555ZUwYcKEPW43b968MHDgwFj8WbhwYVi2bFmYMWNG6NKlS7jlllvCzJkz97Y5AJSB7AcoHtkPUNDiT17Muffee8OFF14Yxo8fX+d2q1evjtvlHnjggTB06NDqdRdccEHsDZS76aab9qY5AJSB7AcoHtkPUNDiT34a1+WXXx7222+/cM8999S77eOPPx62bt0aT/U6/fTTd1t/xRVXxGVVVVVYvnx5U5sEQAuT/QDFI/sBClz8+epXvxoWL14cvvnNb4a+ffvWu+38+fPjcvTo0bWuz/fPTwfbeVsA2h7ZD1A8sh+g/dunKTu98MILcdau448/PnzmM5/Z4/ZLliyJy7znT10GDRoUVqxYEQtKezJlypQwderUBrV10aJFDdoOgPrJfoDiac3s950foBWLP9u3bw+XXXZZvJ4XYDp02HPnoTVr1sRlPrtXXUrr1q5du8f7W7VqVTxFDIDykP0AxdPa2e87P0ArFn9uv/32WHjJB2kePnx4g88TznXu3LnObfIZv3KbNm3a4/316dMnjBw5ssE9fxpynwDUTfYDFE9rZ7/v/ACtVPxZunRpnJErH59n8uTJDd6va9eucZkP+lyXLVu2xGW3bt32eH/5ANGlQaL3pLKyUi8hgL0g+wGKpy1kv+/8AK004POVV14Zq/nf//73w7777tvg/Xr06FGjG2htSutK2wLQNsh+gOKR/QAF7vmzYMGCUFFRET72sY/ttq7UbfOVV14Jhx56aLw+ffr0cNJJJ4UhQ4aEX/3qV2HZsmV13ndpivd8WwDaDtkPUDyyH6DgY/5kWRZWr15d5/odO3ZUry919xw1alR48MEHw7x582rdZ+XKlXHE/9K2ALQtsh+geGQ/QEFP+3rjjTfih0Btl7y4k+vfv3/1baeddlq8bfz48aFTp07x3OE5c+bUOo1jbsSIEfVOCwlA+cl+gOKR/QAFLv40Ve/evasHaM6ni1y8eHH1upkzZ4Zbb701Xm/MYHIAtG2yH6B4ZD9AIqd9NVVe4MnPHX7mmWfCsGHDwjHHHBPWr19fPdbPtddeG3sIAZAO2Q9QPLIfoKA9f0pTOc6dOzfccsst4eijjw5LliwJf/vb38KYMWPCY489Fm6//fZyNQWAMpH9AMUj+wHanoosH5wnYZWVlaGqqirsH94RTqg4s7WbA1AWz2azwrrwRhg5cmTsdVk0sh8ooiJnv9wHiurZBmZ/2Xr+AAAAAFB+ij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhij8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAABImOIPAAAAQMIUfwAAAAASpvgDAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhFVmWZSFhPXv2DGvXrg0dQsfQPezf2s0BKIsNYV3YEbaHHj16hDVr1hTuVZf9QBEVOfvlPlBUGxqY/fuExG3evDku8xdjXXijtZsD0CoZWDSyHyiyIma/3AeKbvMesj/54s8hhxwS/vrXv4auXbvGF2PTpk2hW7du4aijjmrtplFGixYtcuwLrIjHf8WKFTHz8gwsItlPUf/tU+zjX+Ts3zn3Bw4cWMjjz/9x7IutiMd/RQOzP/nTvnZWWVkZqqqqwsiRI8OCBQtauzmUkWNfbI5/sTn+xeXYF5vjX2yOf3E59sXm+NfNgM8AAAAACVP8AQAAAEiY4g8AAABAwhR/AAAAABKm+AMAAACQMMUfAAAAgIQp/gAAAAAkTPEHAAAAIGGKPwAAAAAJU/wBAAAASNg+oUA+/elPh1WrVoU+ffq0dlMoM8e+2Bz/YnP8i8uxLzbHv9gc/+Jy7IvN8a9bRZZlWT3rAQAAAGjHnPYFAAAAkDDFHwAAAICEKf4AAAAAJEzxBwAAACBhhSj+zJkzJ7zvfe8LvXr1Ct26dQtHHnlkmDRpUtiwYUNrN40myscp//Wvfx1uuOGGcMopp4SDDjoodOrUKR7js88+OzzyyCNxm7ps27Yt3HbbbWH48OGhe/fuoWfPnmHs2LFh+vTpjkk79sQTT4SKiop4GTBgQJ3bOf7FIPvTI/vZldxnV7I/PbKfXcn+JsoSd/fdd2cVFRV5FSDr169fNmLEiKxLly7x/4866qjs9ddfb+0m0gSzZs2Kx7B0OeKII7LKysqsZ8+e1bedf/752ebNm3fbd9OmTdkpp5wSt+nYsWP27ne/Oxs0aFD1ftdff71j0g699dZb2WGHHVZ9HPv371/rdo5/Mcj+NMl+dib32ZXsT5PsZ2eyv+mSLv48//zzWYcOHWLxZ8qUKdmOHTvi7StXroyFgvwH4gc+8IHWbiZN8NRTT2UDBw7M7rrrrmz16tU11j388MPVBb7aCjmf+9zn4rp8/z/84Q/Vt8+YMaN6v8cff9xxaWf+/u//Ph67CRMm1Fv8cfzTJ/vTJfvZmdxnZ7I/XbKfncn+pku6+DN+/Pj4I/DSSy/dbd2SJUtiYShfv3DhwlZpH0335ptvZlu3bq1z/Te+8Y14bPOeQNu3b6++/dVXX806d+4c182ePXu3/SZNmhTXjRw50uFpR375y1/GIu+FF16YPfjgg3UWfxz/YpD96ZL9lMh9diX70yX7KZH9eyfZ4s+6deuqe3HMmzev1m3OPPPMuP5LX/pS2dtHy6qqqqo+/Sf/wV8yderUeNvgwYNr3e/Pf/5z9X7Lli1zmNqB/DSuoUOHZvvvv388fvUVfxz/9Mn+YpP9xSD32ZXsLzbZXwyyf+8lO+Dzb3/727Bly5bQpUuXcPzxx9e6zejRo+Ny/vz5ZW4dLW3z5s3V1/NBvktKx7p07HfVt2/fMHDgwBrb0rZ99atfDYsXLw7f/OY34/Grj+OfPtlfbLK/GOQ+u5L9xSb7i0H2771kiz9LliyJy8MPPzzOAlWbQYMGxWX+w5G0/OhHP4rLfDavAw44YLf3xeDBg+vc1/ui/XjhhRfirG15gfczn/nMHrd3/NMn+4tN9qdP7lMb2V9ssj99sr95JFv8WbNmTVzmU3jXpbRu7dq1ZWsXLa+qqirce++98Xo+FfzOvC/SsX379nDZZZfF61OnTg0dOuw5zhz/9DnGxSX70yf3qYvsLy7Znz7Z33w6pN79r3PnznVuk58Sltu0aVPZ2kXLWr16dbjwwgvDtm3b4vKiiy6qsd77Ih233357/MD/x3/8x9jDqyEc//Q5xsUk+4tB7lMX2V9Msr8YZH/zSbb407Vr17jcunVrndvkYwLtOiYM7debb74Zxo0bF15++eVQWVkZHnrood228b5Iw9KlS8NNN90Ux2eaPHlyg/dz/NPnGBeP7C8GuU99ZH/xyP5ikP3NK9niT48ePWp0A61NaV1pW9qv9evXh3PPPTcO+Dds2LDw5JNP1hjrp8T7Ig1XXnll/Cvf97///bDvvvs2eD/HP32OcbHI/uKQ+9RH9heL7C8O2d+89gmJGjJkSFzmvUDyU4BqG/R5+fLlNbalfdq4cWM4//zz40xO+bGcNWtWOOigg2rdNl//q1/9KixbtqzO+/O+aPsWLFgQKioqwsc+9rHd1pVO43zllVfCoYceGq9Pnz49nHTSSY5/Acj+4pD9xSL3qY/sLw7ZXyyyv3kl2/Nn5MiRcbyf/NSu5557rtZtfvnLX8bliSeeWObW0Vzy3h/jx48PTz/9dBgwYED4xS9+Uf2DvzajRo2Ky3nz5tW6fuXKlWHFihU1tqVtyrIsnuu96+Wtt96K63fs2FF9W+n0T8c/fbK/GGR/Mcl96iL7i0H2F5Psb0ZZwt7//vdn+VO89NJLd1u3ZMmSrEOHDnH9Cy+80CrtY+9s3bo1O++88+Ix7NevX/bHP/5xj/u8+uqrWadOneI+s2fP3m39pEmT4roRI0Y4PO3Ugw8+GI9h//79d1vn+BeD7E+b7GdXcp+c7E+b7GdXsr/xki7+PPfcc1lFRUW8TJkyJduxY0e8/S9/+UtWWVkZfyBOmDChtZtJE7z99tvZhz70oXgMDz300FjMa6jPfvazcb+BAwdmf/jDH6pvf/zxx7MuXbrEdT/96U8dlwQ/CHKOf/pkf7pkP7WR++Rkf7pkP7WR/Y2XdPEn953vfCcWf/Ifg4cddljs0VH6gT906NDstddea+0m0gQ//OEP4zHMLwMGDMhOPvnkOi9VVVU19t24cWN24oknxn07duyYDR8+PBs0aFD1/V177bWOScIfBI5/Mcj+NMl+aiP3KZH9aZL91Eb2N17yxZ/crFmzsnHjxmU9e/aMhZ8hQ4ZkN954Y7Zu3brWbhp7+Y+9IZc5c+bstv+WLVuyW265JTv22GOzbt26ZQceeGA2ZsyY7LHHHnNMEv8gyDn+xSD70yP7qe99IffJyf70yH7qe1/I/oaryP/TnGMIAQAAANB2JDvbFwAAAACKPwAAAABJ0/MHAAAAIGGKPwAAAAAJU/wBAAAASJjiDwAAAEDCFH8AAAAAEqb4AwAAAJAwxR8AAACAhCn+AAAAACRM8QcAAAAgYYo/AAAAAAlT/AEAAAAI6fp/lnjksgvLUvEAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(ncols=3, figsize=(14, 10))\n", "\n", "im0 = axes[0].imshow(mask)\n", "axes[0].set_title(r\"$\\mathtt{n\\_erosion\\_threshold} = 0$\", fontsize=14)\n", "\n", "n_erosion_threshold = 5\n", "selem = np.ones((n_erosion_threshold, n_erosion_threshold))\n", "mask_er = binary_erosion(mask, selem).astype(np.int64)\n", "\n", "im1 = axes[1].imshow(mask_er)\n", "axes[1].set_title(r\"$\\mathtt{n\\_erosion\\_threshold} = 5$\", fontsize=14)\n", "\n", "n_erosion_threshold = 10\n", "selem = np.ones((n_erosion_threshold, n_erosion_threshold))\n", "mask_er_more = binary_erosion(mask, selem).astype(np.int64)\n", "\n", "im2 = axes[2].imshow(mask_er_more)\n", "axes[2].set_title(r\"$\\mathtt{n\\_erosion\\_threshold} = 10$\", fontsize=14)" ] }, { "cell_type": "markdown", "id": "dd64d324", "metadata": {}, "source": [ "This means by using increasing values of *n_erosion_threshold* for a feature detection we will get lower values of *num*, which will match the number of **True**-values in the masks above:" ] }, { "cell_type": "code", "execution_count": 6, "id": "2e7bdf61", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.750646Z", "iopub.status.busy": "2026-02-02T20:12:30.750530Z", "iopub.status.idle": "2026-02-02T20:12:30.830765Z", "shell.execute_reply": "2026-02-02T20:12:30.830222Z" } }, "outputs": [], "source": [ "%%capture\n", "\n", "date = np.datetime64(\"2022-04-01T00:00\")\n", "input_data = xr.DataArray(\n", " data=gaussian_data, coords={\"time\": np.expand_dims(date, axis=0), \"y\": y, \"x\": x}\n", ")\n", "dxy = input_data[\"x\"][1] - input_data[\"x\"][0]\n", "\n", "features = tobac.feature_detection_multithreshold(\n", " input_data, dxy, threshold, n_erosion_threshold=0\n", ")\n", "features_eroded = tobac.feature_detection_multithreshold(\n", " input_data, dxy, threshold, n_erosion_threshold=5\n", ")\n", "features_eroded_more = tobac.feature_detection_multithreshold(\n", " input_data, dxy, threshold, n_erosion_threshold=10\n", ")" ] }, { "cell_type": "code", "execution_count": 7, "id": "ddefe3e3", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.832220Z", "iopub.status.busy": "2026-02-02T20:12:30.832086Z", "iopub.status.idle": "2026-02-02T20:12:30.834415Z", "shell.execute_reply": "2026-02-02T20:12:30.834043Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(332)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features[\"num\"][0]" ] }, { "cell_type": "code", "execution_count": 8, "id": "c8c0733d", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.835625Z", "iopub.status.busy": "2026-02-02T20:12:30.835545Z", "iopub.status.idle": "2026-02-02T20:12:30.837640Z", "shell.execute_reply": "2026-02-02T20:12:30.837210Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(332)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mask.sum()" ] }, { "cell_type": "code", "execution_count": 9, "id": "5909a742", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.838714Z", "iopub.status.busy": "2026-02-02T20:12:30.838640Z", "iopub.status.idle": "2026-02-02T20:12:30.840681Z", "shell.execute_reply": "2026-02-02T20:12:30.840314Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(188)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features_eroded[\"num\"][0]" ] }, { "cell_type": "code", "execution_count": 10, "id": "33683d04", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.841845Z", "iopub.status.busy": "2026-02-02T20:12:30.841769Z", "iopub.status.idle": "2026-02-02T20:12:30.843811Z", "shell.execute_reply": "2026-02-02T20:12:30.843448Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(188)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mask_er.sum()" ] }, { "cell_type": "code", "execution_count": 11, "id": "2cb7802e", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.844926Z", "iopub.status.busy": "2026-02-02T20:12:30.844836Z", "iopub.status.idle": "2026-02-02T20:12:30.846881Z", "shell.execute_reply": "2026-02-02T20:12:30.846529Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(57)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features_eroded_more[\"num\"][0]" ] }, { "cell_type": "code", "execution_count": 12, "id": "98b07f0f", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.848023Z", "iopub.status.busy": "2026-02-02T20:12:30.847949Z", "iopub.status.idle": "2026-02-02T20:12:30.849997Z", "shell.execute_reply": "2026-02-02T20:12:30.849600Z" } }, "outputs": [ { "data": { "text/plain": [ "np.int64(57)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mask_er_more.sum()" ] }, { "cell_type": "markdown", "id": "a953d755", "metadata": {}, "source": [ "This can be used to simplify the geometry of complex features." ] }, { "cell_type": "markdown", "id": "b16cfa52", "metadata": {}, "source": [ "## Minimum Object Size Parameter `n_min_threshold`" ] }, { "cell_type": "markdown", "id": "efda10e4", "metadata": {}, "source": [ "With *n_min_threshold* parameter we can exclude smaller features by setting a minimum of datapoints that have to exceed the threshold for one feature. If we again detect the three blobs and check their *num* value at frame 50, we can see that one of them contains fewer pixels." ] }, { "cell_type": "code", "execution_count": 13, "id": "1ca6e991", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.851208Z", "iopub.status.busy": "2026-02-02T20:12:30.851118Z", "iopub.status.idle": "2026-02-02T20:12:30.917030Z", "shell.execute_reply": "2026-02-02T20:12:30.916712Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAKrCAYAAAAebZGTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXxJJREFUeJzt3QuYHFWZ//FT3T0zuZCQIMQEAgTCzQBCErmGyC1Bs2BQCMrKI1cxyKNGlucBHhDB/4Ig8riwuCuJIEFYBAO4wLqL3IW4sBES4kazBEICCAGBBMhlbt1d/+ec0KfeU1Onu7qn59rfz/M0qamuqq7uGeZM/eo95wRhGIYKAAClVIZPAQBQQqMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAAAYuI3Ck08+qU444QS1ww47qKFDh6p99tlHXX755Wrz5s19fWoAMOAFA2nso5tuuknNmzdP6VMeP368aRj+8pe/qPb2dvWpT31KLV68WG233XZ9fZoAMGANmCuFF154QX33u981y/Pnz1evv/66Wrp0qXr11VfV1KlT1cqVK9W5557b16cJAAPagLlS+OIXv6geeOABdfrpp6vbb7/dee7ll182MVKxWFTLly9Xn/70p/vsPAFgIBsQVwqbNm1SDz/8sFn+xje+0eX5PffcUx1zzDFmedGiRb1+fgAwWOTUALBs2TJz36ClpUUdfPDBidtMnz5dPfbYY+q5556ry2tOmDBB/e1vf1NDhgxRu+22W12OCQB9Yc2aNaqtrU2NGTNGrV27duA3CqtWrTL/7rLLLqqpqSlxm4kTJ5p/X3rpJe9x9L2IBQsWpHpNfc9CJ2utra1qw4YNNZ03APQn+g/dSgZEo7B+/Xrzb7nKotJz5X6Br1u3ztycrkZGZdVwNaLMFkHlgwQp9xGrA+cL5wnfDp7tk18rjK/3vYZnvXMjyrd9T/HcBgvkamdZfOHsGnr2DVPsGz8Pz3Zim7DLAcqfk2/zMjsDiTarjaqoCib5GBSNgr7s0Zqbm73b6GhJ03/Z+4wbN05NmTIl1WvqaiZ9LN0gHBLM8P+yC5JvywSZwL+NeC6Qx81mo/ViWa5X2ehYQS6XvE1TtD7MJR8nFNsYuei4YVO0XVGuz0bnGvrWy2XxPkOn4XBfOtpIbFLmF3BQiJYzYjnIy22KYr1Y7hTLhYJYHy2rvFyfT1y/9evouVC8nlyvxGuEYtlpLOT6omwUxDGFMMU2XXei8Whk/xM+pjaqD1JF4QOiUSi1bh0dHd5t9D0HTXdo85k7d655pKHLXO1VRcq/fr0NgVxfz4bAWfb8wpcNhPhlL5e1onxO/MIvNsllzy9/sVzMyUZBJS6nEW8UAvG7z2kUnIYg2iYjfvlnZKOQE+s75fdIfMZiWV5ROd+3rSvEomggfO9JLMuGQH7fQyUbiEzlX/7y5yxtAwEM9Oqj0aNHOzFSktJzpW0BAIO0Udhrr73szd/Ozs7EbVavXu1sCwCo3oCIj/R9AH0/QUdES5YsUdOmTeuyzTPPPGP+Peyww3r+hMQluxMZ+TaPxw4innCXg/pHRs3J8VGxORYfNUevV2iuHB8Vxcs5kZE4rHNPQaYcvtszvhvFsfhI3kfIiPheRkmZziAxJsqI+CjMiihJRGDy+9XleyfPybde3i/w5GbeKEm8nrOriJKCTDH5/kI58n1wfwED/Uphm222UZ/73OfMclJJqe7R/MQTT5jlOXPm9Pr5AcBgMSAaBU2PhKr/irrjjjtMw1AanUOXmf793/+9GeJCD4VxwAEH9PWpAsCANSDiI+2ggw5SP/nJT9Q//MM/mAqiq666Sm2//fZ2lNS9995b/fznP++5E/CUnnq3kVGQjIjMdiKeEDGRW3HkWfZFRs1NyTHRkFxihVGhJeOPj1qSI6OC6Dfoj5KUpzy1ypLUWCFNIGKSoCCjoWibTD5an+2Mts92iFipQ2zjqaDKiO9P/FvnE/jKZ+VqlU+uFHKyNbFvUcREafojxH9GqUbCYL5S0PQoqY8++qiaNWuWmT9BNwi77rqruvTSS9Xzzz9vGgkAQANcKZQce+yx5gEAqL8B1yj0B2k6qTlVK7EKFqeTmrO/6MgmcwsZGckeyp6OaU5k1CIjI7Es4iLz9RAREzVHy/kWERM1e+KjJl98pCr2aJbry1UfZURkJPqJufGRiIyKIiaS55T1VEplU1SRZcpVGXnW+4fMkDGROI78WZE/A3JfKpHQgwZUfAQA6Fk0CgAAi/iooiBd5VGcr4NauU5qshIpzWB3ns5o3shIVBXlh7rnJGOiwtaxBT9eTl7vRkme+CgXJkdJKpkT4MSrj5xxjTzVR2J9sT35nIpO57qMp9zJ7djnJTvUyShJVA0FvmoiERWqrCdi8nVqSzsyKuMioQZcKQAALBoFAIBFfJSSd4yjFBVHXcbPqbaTmtNhLZs4rpGMjwpivCJfZJQX1UZbvxbLQ0VkNCQ5Miq0hMnxUXPoiY9kzzSVovNarGJLdEwLnMgoWs62J1dEhU7FkTgn53vqDDQklkVFWGzMoFBkYr7IyJkfQXweTuWSeK8yQpTHrLYSaev+aTq8MSYSXFwpAAAsGgUAgEWjAACwuKdQC1+Japp5Esz+QeWey7I81VeGKuc6kKWnQ7KJ9w6c5dispfI+gnxO3lMoDBH3EeQ9hRaRZTeJbF0sZ0SW79xiCZKzdblsXkNOnSmWnZ7L4vPItvvmdQhS9KZOnuIykPWsmpwfWk5vKtYrMX+DLEl1amPl/M7yPoD82aA8Fb2EKwUAgEWjAACwiI+6KVUZajxukqWnvshJbBPmMhV7McsB7txeyDI+So6Ltn4tlodFEUZhqIiJxHI4JOptm2mJlpuao9Hqmpqi9TlRLpmJlU7aY4o8Jy9KLbXOzuhHtbNTvO92+Rkkx2/uXA7J9bBufCS+jyLGkvM1mOfE68lB7eR0oU5klJfdspN7sjvlqbIXs29aUBlNisOb/WVUlnbaTjQ8rhQAABaNAgDAIj5KwVyG+6ba9PEMemeO5xkz36k4cno0J29TFD11i82V50NwYiURJXWJjJxlUT0zNMonmodGXYmHtETLw1o6ovW5KEpqyeYToySpKGKbzlilT2tn1EV5S0e03NoUdbPuEIMIFgPxoy0Hn1OVY6JAxD8yMsoUYt9HmQYVkqujws7k76kbK8mfBxktFipO2SmrkpiOE/XAlQIAwKJRAABYxEfd5aseknFTueojGRk5UVKQHBnJDmui4sidHjN5Ck23I5p7Ss5zMjIaJiKgYVFMtM3QqHfYyCFtdnnb5la7PCwXbT88F23f5MxdEOkUEUmHnARBKbWxM5rMYWNTdLIfNUXntzEbrW8ViYzzaiKWcuMfOd2nZ74GMQ2o+VrGd3JeDBH3ZcWykhGTjBDlz4OMlcTPTSA7+YlzqGmeBR8GxwNXCgAAifgIAGARH9Wgy/wIFdaXG/vI7bAm4yNf9VFyTFSQy868B571Yhwjc6whIrYYklxlNGJYFBNtN3RLtNwilpuj5VFN0fKwTFSV1CRymKyIPDrF/ARbZO6llPqwaWjicktuuF3OOD3QIq1y+gFRZVTMi45vnjkanGlHxfpykV1GjMckv3eBrEqScZOcK0LGQb6OaWl/zmKd2aL9k8d2AsyPER8DAKCERgEAYBEfdVeKS/kusZIvInCipMrVR06UJCth5PSYzlSZ0XJBTJtpjiWGv24aEsU7w0WV0eghUWXR9kM22+VPtnxkl8c0b4y2z0XbDM9Exxki59MUZHy0uRhVG2kfZofZ5fcy29jljIifZOe3ghg7KZ+PltvluEkdmcTIqNAZHScrI6LY/y1FZ7rR5MjIGao76xkyXf4MeMbSKjfGkZesXhIdBhkHCeVwpQAAsGgUAAAW8VFa1Y535IsBymznRBCe2EEOgeNERrkU60VkFMbio0AMf93cEsVHI8RYRqOHbEmMjMa2fBitz0XLn8htio6TaetWfCTjpybZ68zZX8REheiNt4tht/NiWcZHRTEEd5jzxG9d4qM03ztf50ZfTJSiA2QmxThIQI24UgAAWDQKAACL+KiSIGHsIt94R2kF1VWbeCeelymCGD5Hxh/OsrONGx9lm6LqlJYmMRR2UxQfbdsUVR+NFh3TdshFFUdjxPIO2Wh5hCjvaRIVQyJ1UZ3ilDbGIqZmT2TUJt7gpkIUOW1qEsvNUdlVa1OUBxXEUNZFUTHkxG/yM4uNwO39/OWPhPx+yZn5fFEj0Me4UgAAWDQKAACL+Kie0lSLpIyPnE5tQeWYwrssR2WWkVHWjY9yuSieGSKGox6Wi+KjbbJRBdC22ShKGpUVYx9lo4qjUWK8oxHi/bQEyTOhdYpxeJo8FUpaWzaKgD4qijGRcm2Jw3bLGeByueg1OmTnM6cjWnKntKLMurp85p5KJLlLpvYIUXmHywbqiysFAIBFowAAsIiPBkCnOCeCkFFSkCZKCpP/BIjFRxmxXZMYJ6dZTB4/NBtFMi2immhYEMVKwwMRPYnzGyYio5Ygin8y4g3lZYVRMYqetDYx3PawsD15TCVxTvJc5XvIivemPJ+N77OUn32X7YIqo6FqK45SjH0UH2Or2zOxoSFxpQAAsGgUAAAW8VE/EqbpyOZs4zmQL2JyIo7Y2Efia7mck9FLIGIlERPJjmVyKOsmca5NIj6Sy1JRjOPTFOsw2CQyk6wqJi7LMZEy4lzl+5Gzs8n1vojOUS7x6Ub/M1m55J29D+glXCkAACwaBQCARXzUX8khkYsi8pDrfZlFmJwSyeVydSmhyFKKosSm4Fv2nIecEr4gzrsg45/4uFL2dd0zTDO9vDwned7ucSt/ZjUV7dSr0Cf2voHexpUCAMCiUQAAWDQKAACLewp9zXu/IAXP/QI3Hw889yncbL1YjP4+yMtlkc3L6TLlclsY9VDeIuY3GB5GvZKbVL7ie+gUZadtYnC8rV8nv55cbhev7Zy3mBRB3lMI5WeQ4j5M/MZGrKpXHCusajkQPwNeKbYJuR+BOuBKAQBg0SgAACzio3qKRR5WMbbed5nvxAvRouicq+SYcc76Yort5UBqeTc+KuSjvw86CiKqyUfxzOZ8NMXlxqYh0XIhmtNgeCAiI3kiYkC7NhErZUUP3g7x/jeKCEv7oCheW8yhIF97SyGadnNzPlpuz4tYKR+9t7AQvUamEFT8zOIxj/za+cxTxk+1IiZCT+JKAQBg0SgAACzio0rCUiwk57WUWYGIYbJygH2ZFcQGgEtTkVIQA7p5KlXSxETOsoyMRFwSj4/aO3OJMcxHIj7a0DncLg8T027KQemkzrAtMVbKipysU1QMfRRGr6V9UBhml/+WH2mX1+ej8/igM9pmU2dLcnzUKb4XnSIyEsVRgS9KilcfOZGd/N55vqfyey2+v96fhzRVSYBWaSDFKn6UuFIAAFg0CgAAi/ionuTlvi8SKLedr4JFbi8LekRMISMPUejjLkczVKpMh/v3QLEjilXaO6Ifi43tUQwzNBfFMy1inoUmsSwHopMdy0ZkWj3zL0RvqFPEbJtFtZG2Pr+NXX4vP8Iuv9sRLW/oiCqRNnVE+29pj86j0BmdXyCWMyJKcj4nz2cZj5YyeV/EV11HNiealBGkr7INg0PQf+bR4EoBAGDRKAAAej8+0h1unn32WfXggw+qxYsXq5UrV6qPPvpIjRo1Sk2ePFmdccYZ6qtf/WridIQTJkxQr732Wtnjt7a2qiFDog5VdSdjgGzQrc5GTjWRrEIR0UEolgMRTcjISEYWzrKMP0Qsku0Q8yS0u+dYbBFjBbVHPxabclEMk8sUE6fdlNqaxThIxebECqUWeYKCHE9J7qt9mI+iq/Wi8um99mj5/bZoeWObqD4S8VHYFr1Gtl1ERuKzcT8/+Rm75+tERk6Vl6cSyVOVVDZqTCIjpnjHSElETiGVTH0n6D/RUL9qFJ544gk1Y8YM+/Xuu++udtttN7VmzRr16KOPmsevfvUrdd9996mWFjdPLtlvv/3Utttum/hcJsNFDwAMqCsF3Qh897vfVaeeeqoaM2aMfe6OO+5Q5557rvrtb3+rrrjiCnXttdcmHuOmm25SRx11VG+dMgA0nF5rFA4++GD10ksvqaam6FK+5Gtf+5p644031GWXXaZ+/vOfqx/+8If9+i9/GQfJKMhZ71QMxSIBX1wgO6zlxbJYn+ksJkdDIvIoiOVslNoomciINGfr161iKstM9GPRJr4NH4n4SJLDVLcWmhLHRxoqXlDGUFlZfSTio1YxjpG2WXy9oT2Kkj7siF7jw9ZoeXNrtH2hNXo/QbuoOBLxUbY9+bPxRUnmaxkTye+X/H4765O/p/L7HnanIxsR0eCMiYJ6/C5Mf2699pt35MiRiQ1CyaxZs8y/69evV++++25vnRYAoD/2U2hri4ZBGDo0qjeXbr75ZnX99debm8pjx45V06dPV6eddpoaMSKqVQcADIJGQd9k1g444ABzVZHknnvucb6+66671OWXX27+nTlzZsXXmD9/vlqwYEGq89HVUd1SruORr3pELsuOaTKCcCILWXEUJkdG4uJMVtsUYxVUIrlRYTb6QowurbaoIYkztcmhtjd3RrHNB01R4z5E5FsyPpLkrGhtIobSWsUQ3ps6xBDZ7dFym4iMOluj7YNWUXG0JTrvbJusPopeK9uR/Flm4/GRU5kkvi+dhcTIyBnvyOmk5omJxLIvVqppGG06wvVdTBTUL5wJMlWcR/KQZP23UVi6dKm5CtAuueSSLs8ffvjh5n6DvjLYZZddVEdHhylr/f73v6+WLVumZs+erf7whz+oKVOmlH2ddevWmdcCAPTTRuGdd95RX/rSl1RnZ6f5V1cmxekrAWnYsGHqhBNOUMcee6w64ogjzC/6iy++2JS1ljNu3LiKDYe8UtAxFQA0kiDsw2mcPvzwQ3X00Uebv/anTp1q+jL4oiOfhx9+2Nyk1tVK77//vukMVw/6fHRjM0KNVodkj3Mv1US84nS2k+vFsmqKtb05UQ0jb76Ljl9hS/JyYahcjo5TGCo6nw1LXhb9v1R+mHvp6T4X/UgUhkbLxaEiFhHLTS1Rr67m5mh5SFO0nMtG22c9M97LtZ0ikooP590hl0VHu6LomOaLjHKt0fvObVGe5ehMmsRyrtWNvXKbo/eUaxXjObVFUVmmTfR4a4+yqKBdlDV1RsthXmwvlsNOsb4QvZbzv2+XKjdP5zVffNR3vwpUo8dEQTVRUA2v8Vznf6mN4QbzR/ELL7xQdts+q/vctGmT+vznP28ahH333Vf97ne/q7pBKEVLWrFYVKtXr+6BMwWAxtEnjcKWLVvU8ccfr5577jm11157qccee0x94hOfqOlYzc3RzcW8/CsLAND/7yno0tMTTzxRPf3002ZMo8cff9yUl9ZqxYoVdnmnnXZSvUJWhYhmNVVHttj+MgpQRRGZ5MV6EUVlRGVLKCqIwpwc40icR1aeq9g+drXqfi2+cCaej86jKIbq7hDDcHc2i45vOdHRLiuWM574SCYhsuwpNjNcUc6eJl4705ZcWZQVkVE2qnxWWXHLKNsmYiKxnG1PXjZfi46ETsWRqBaT30fvOFf16rAWi4UY76gXIqMgU59oqJaqpKoip37YeU3TN5NPPvlkc2Uwfvx4cw9B/9sd1113nfl30qRJ3T4WADS6XmsUCoWC6Wj2n//5n+bKQDcIeiykSnRnNT3mkb6JLOmv586dqxYtWmS+/sEPftBj5w4AjaLX4qNf//rX9he4HuL6rLPO8m6rGwE9nLb217/+Vd14441q3rx5Jm7aYYcdTKmoLhnV9xB01dE111yj5syZo/qc7JAkL0O7dF4TEVBRTh4vogY59pOMmOQMYTmxTXshOSYSUY0TEXW5mgySZ33zzOhWEDO1yWG3w6Zo53xODjcuD+qJReQJxgtkxHnIGdPc4a/Fy7V7IiOxnGtNExmJcZrEZ2zOo0PGR8lRkvO9y3vWe5bpsNZPK4yCTP1jojL7Jk0nUFbSuHH6EGE/axTa26MRx9auXWse5UpVS3S/Bf0/x5IlS9Trr7+uli9frrLZrBl6+8gjj1Tnn3++OvDAA3v8/AGgEfRao3DmmWeaR7UOPfRQ8wAANECP5oFCVnIEcuweMQaQMyObrESSVSfmcjBTuRJJVt+IyEheSsqISa7PZkQEEaS8hSTTHZmCydhGRklyfCVnrKXAM56SJ8ZSnlRJxGrm9ZyZzVTFGdPkUNjOsqfKyF2WkZGomhJx0dbXKyRHRjIGlJFRtRVHaWZY6+44Ro3SYa2OMVGqyChFTOSNhcpNG1BzlNRPq48AAP0bjQIAwCI+qijceome4lLS7bDmqUSKx0SyUshXiST2D2VkJJadK84gzTfYfT8yMsrImEgM4S1nepOzuBVFhCP734XixWVFlPdK1hNhma9lfFTwzYzmGf7aiZLksNieyEgsZ0TFUaYjXn0U5ViBGJvIiYx6ouLI06mNzmp9UFkU1CkmKneunpjYv33QrY+CKwUAgEWjAACwiI96uRJp65difxkjyGs8z7Jz+ei5Jsx41stvdrz/WCCjK7GcySdX9zgVR+LAxZyn+sgTbznK9G8LfNVHzuxzybOkObOniQqijKeTmoyJZIc1GReZr53qI/Gcb5hrWYXmjIPUQxVHjThEdh3HLgrSREbdiYl8sVC5TnCe1/ZvH1T9uXClAACwaBQAABaNAgDA4p5CWjKf7U55qlnhuS8gcuZQ1GQGTqCeIsMUMk6JozhMrKwxEPMVZMRyoVmsFzm97LlcFPdPimIQPOc+glhO1aO5zD0Fp3ezLJnNJ99TyIjS04wYuM69v5BceurcN4jdU/DdO3Cm0XR6NHvKUOW+lKH2vLT3EYIqy027c++g3L0C5/dFdfcUotfjngIAoAbERwAAi/ioF8pTnRglPm2niA5SVGp6t0lzcShPL4xNd6kKWU8v5jCx3NSJj2QZqoyJ5HShZedyKG1UrmRWLMvzE8veUlUZE4n4yDegnW/ZiYviZagiJgqdstLqei67JanJU216ey6XK1Vt9DLUtD2Vg25ERmLa3FQxUTaTKhbylqGn7RGd5nmBKwUAgEWjAACwiI96qhIp9PR07jL6XKZyT2c5aF4+RZTkiSZkFVQoR66LVSMVRZQkp5nMNGUSB8ELxbSgTnzkGQQvHqdF5yC+iMdHzmeTfN6B07tZVG+Jqi5nfZr5EHwVRrHtQlFxJKuP/D2aU8yzkEZ351No9Mgovk0m+biBiIZSVRbJ7eUxM8nrnel34/MpOM9VHvXAPXF6NAMAuoH4CABgER+loS/pfZdnTlWInOsgqK1jW2kbWckQjy2SoiRfZCQrnYrJkYohqniy+WJyNCQ7uMlIS1RRyIojJzJyLnVVZV2qj2R8FHrWy/cnDpD3REaFFFNoOstu9ZFTZZQmMvIMduedKyFNxVE5g63iqKciI9+hsikiI18FkVwv/1/2VCs5x5f7xl7PjWRTLJe2TzN4XmmX1FsCAAY9GgUAgEV8VMuleIrLWG+ntjId27xRkvKMWZSV0zWKS1rPeTsRUzw+kq/tVB9lPFVGnrkj5NShTue16sa5d87VnF9yZCRjL+fzk+udqS9FlOSdKjM5FnKioLJjGXkiI1/ntXp2UmtE1U6hKcU7imWrjIycKqNsim08EZOMiMrER26Ht+T/19w3FJQdbywJVwoAAItGAQBgER/1Qqe2eAyQZoykblUlKVEJ46t0kpe08edkpyxxuSqjpMB36SuH9nam4EzR0aZctYzTec0zbpAntgm8cY4nJvJEPl0iNxn1VNsxrScio8FWbVTL9Jre46ScQjOoU2SUE79aPfs6MVHOEzHFoyEnZqpcfWRjI8Y+AgDUgvgIAGARH/VUJVKZ8ZGqHW7bGyXJqMb32rLzmjzX2DEDOZS2jKjynk41ng413mF+q40B4lGIdzynFENQy2jHiZt8FUMy/knefuuhwrpERl6NGhmlUW0ntTRDX9cSGeWyyduI9WGabWSlXizaVbLqT27n6cjmrLfnTOc1AEANiI8AABbxUW9HSWUqk1LN3CZnapOVDM7k7xlPZ7esN3bwRUvOpbYcbltGV75ZpupZReKNYSrHR+5MaGHlmMhXVRSP8XyvR5VR91Tb0bGqCez9s6V12S7IVNcxzRcZNeWSoyEZC4nlYlMsPpLjjDmzHHqW5cfx8frESMmDKwUAgEWjAACwiI/6YnykFJ3cfFGScxgloiT5er6KI1m1Ez+/TIr9ZZRSrgNQ9ETy+njFR4mnyqoLeX5pqnvSxDm+mMhX0VTutZ2N6JhWN2mGvPZUHLnblPnZzcpKoUzljmlpIiMRB8n4KBQzGTrbiIhIK8poyRMfyT/v3Vjp438Z+wgAUAviIwCARXw00KIkKU2FUppYqVzFUooZq7yzQVUbA6QdEtrXYcsXH4n1buRTTHFMT1VRudeTp8pYRv2v4shZLhOjZj2dzuRMg2kiI8+yjIVklFQUy1u3Sx4i242SVPLvgtJnUMWf/1wpAAAsGgUAgEV81BvisUOV4yX5IohUFUopxlAqX3HkGRfGu02dhjoux/N5dIl3qqwmSjyOL/4xGxIZDfiKo/jxM6KTmi9y8sRKTmVR1hMZOcsiMmpOrjDaup2MjFRylJT1xbkfr2PobABALYiPAAAW8dFAGS+pHhVKylNVVK5iyXPeTqjii5LSdFjzqaEjW5pj+SOm6mc881YWpdxfnFTlbZBK1RVHslNal+eyidt5Z0zzVRM54xolR0aFluSqonLxUdEZbluVrURyKpIq4EoBAGDRKAAALBoFAIDFPYWB0vO5HmWrkqeEtey9hzT3HXzkfA/dHRDPo+r7Bd0pL025v3977iN0qxdzvcpQYyWpgfzZ9N2f8E2j6dxTSB74znsfQdw3KIjlrs+p5PsIznLXuRUYEA8AUBPiIwCARXzUn/gihTqVrbqbpIw7vHM5eM41TY9mGSXVIm2kU0W0U7fy0i77EBP1Kl806cw9UGb62Eymchmq2Eb2XC565kRwylDl+lxyZFRojsdHnmVnoDyVvFzahJJUAEAtiI8AABbx0WCvUKoyVuq6S5WRSaFOg5zVEtV4D1X/uKn8/kRG3ZLm50NuXu20r+Wm48z4ekHLyEjsIyqOZBTlm0LT6Z3sXXZPScZJ8rlQ9m7O+eZT+Hgd03ECAGpBfAQAsIiPBnuFUi2xSA0xUzV8kVRPvFbCi9TxWMRE/U61Hda6xEdB8rwE3uoj2XktSDGFpkpclp3S4vGRt/rI15FN/qlfWiY+AgDUgvgIANA3jcKVV15pLunKPW6++ebEfTs7O9WPf/xjdcABB6jhw4er7bbbTh1zzDHq/vvv78230H/pKCPpUdOxitU9qj58mPjoj+fq/VyJjnqFHu+o9DCxT+mhvy490pDbO4+M+wiCig8dB9lHpvJDR0mlh9y3KB46/ik95PZmn5yyDx0ZlR6+9cWmsMujmvioT+4pjBkzRu25556Jz40bN67Lura2NjVz5ky1ePFilc1m1b777qs2b96snnzySfO4+OKL1bXXXtsLZw4Ag1ufNAqzZs1SCxcuTL29/qWvG4TddttN/dd//Zfae++9zfoHH3xQffnLX1Y/+tGP1LRp09QXvvCFHjxrABj8+v09hXfeecdGSrfeeqttELTZs2eriy66yEZTqDL+qFcUUm2EU89Ht867hz8XdCWjmP4gRVykYpGRdxsRGZnfrB8/QvEoZqNHmPKhq5RKD+/6XBg9dLT08cO+1mCqPtJXAx0dHWqPPfZQRx99dJfn586da/5dunSpWr16dR+cIQAMHn0SHy1fvlx99atfVW+//bYaMWKE+vSnP61OPfVUc68g7rnnnjP/Tp8+PfFYO+20k4mV1qxZY7adOHFij58/AAxWfdIovPjii+YhrwauvvpqNW/ePHX99debm8klq1atMv/qKwUf3RDoRuGll14q+7rz589XCxYsSHWOK1euVA2l2qikv1z+S8Q9jSnNcNm1CCrP6CY7uMmIxlmfSbFedoITHdG67i+Xw4qd16qJjfqkURg7dqy5B3DSSSeZX+T6KkH/0v/Xf/1Xc9/ghhtuUM3NzebGccn69evNv7oE1af03IYNG8q+/rp160zMBADoB43Ceeed12Xd/vvvr372s5+ZCEhXGf3kJz9R3/zmN9WECRNsOaqmGwuflpYW829ra2vZ19flrlOmTEl9pVDpeAAw2PSbsY8uvPBCdeONN6q33npLPfTQQ+rb3/62WT9kyBDzr77Z7NPe3m7+HTp0aNnX0DelSzemK5k6dSpXFeUQ1aA/qmW4bM92SsZPMjLy7e6JhpxjOhGTWJbr48eXXyeNa2RmQvTsX4qYgnDgVR/p+wiHHHKIcx9BGz16tBMjJSk9V9oWAFCbftMoyIgon8/bdXvttZf595VXXvHuVypFLW0LABgEjcKKFSvMv+PHj7frDj30UPOv7tGc5M033zSVR3JbAEgl3gmtWrKTWhA9qhb4H2GKh7tPaB+1nE+/aRR++9vfqj//+c9m+bjjjrPrTzzxRNXU1KRefvllM85RUpmpNnny5LJlqwCAftQo6F/4+iav7rgmFYtF9atf/cp0ZtOOP/54ddBBB9nnP/nJT9qbw+ecc47TF0HfkL7uuuvM8hVXXNFL7wQABq9eqz7SQ1/rjmP6ofsV7LrrriqXy5l7BaX+BbrX8p133tllX/2L/4UXXlDPPvus6fW83377qU2bNtl7CbpySV9RAECvVtEVo8UqCnwSzsOzHDuub9l5GyIvCmp4f712paD7HVx11VXmSmDUqFGmMdC9mvXNZT1q6h133GHiIf1cnC41feqpp0yntkmTJpnqpPfee08deeSR6t577zW9oAEAA+hKQf+yv+yyy2reXzceujd0aVRUAMAg7rwGAN1WFHlOVvbuEuMEieUuhTmhiFvkbIByH18iI7YPip5jymjHEz11Ob78Wh5X7i/HS5IvV/z49aooQeo31UcAgL5HowAAsIiPAAyiyCibHP/EhqNOJfTkRDJKkpGRjICcKCmovL4g1hfcqMeJmZxlsZ3Yv7u4UgAAWDQKAACL+AhA75GRTH+YvS8eEYWeiqNCcvWSL0pyKoPEcqYgNhHLgWfZ7JNPHm5bjuxdlHVUMqLKVN+xjisFAIBFowAAsGgUAAAW9xQA9Fuhk4+LcL4o/p7Nprg34eT9crmY6p6C8txfcEtMk5cz+Wi5mPPdXwgStzcvLW4euPcRxLHkqcqNSh8T9xQAALUgPgIAWMRHAAYP38B3oQhbnOXY38XF5AgoLMjoqpgcJYnYR0ZAYVauj6KdMJs8dl98lD7fWHayzFTGUs5bKu1LfAQAqAXxEQDAIj4C0L/IeKfUJbfc5r75EaqcW8HsX0yOiXyRUZgX2+Si18iI9aGojgpFmZCMeWREFHbp6Z086F7RGVBPJcZSpePSoxkAUBPiIwCARXwEoLHmVnAymNC/f6GY2HEudKKkKLfJdEZxTlF2OJP7ivUy5uk6L6g8XZktqeQYS0ZGMq4qLcb66JXDlQIAwKJRAABYxEcABsTcClWPg1RtRzbzZSC284yRlI8ioyAjXltUHAViYCP51jLO+/T8TR5LtNxpPsUpyZhIxkcJFU5UHwEAakJ8BACwiI8ADBpOTCQjHxnbyGU5ppF5LiOek3NkilhJVhPlPduI5Ywcytobk2W803HKyEg+l5HFVU58JF+v6zEq4UoBAGDRKAAALOIjAAN7HKSwTpVIZv+COKyIg5xObTLDkbFSwoxnWmfyamc2Nyf2ct9nIM5DDr3tzOImx1GS07B9HFc5lVQVcKUAALBoFAAAFvERgAHekS2oTyWS+TqsWH3kVhmJuEkexnfi4jQycmY3p4NabDjvQvS3ezEnXltESfLPe2d8pVKURPURAKAWxEcAAIv4CMDA0MOVSOawhVjPsdJ2ovdXoPKJ28jIyBslybcgz0NERkXZE02fej56LiPiIzcmShgu2+xQqj5SqXGlAACwaBQAABbxEYDGrUSKb5fNJndsk53axPay+ihVlCQri+R41+KcMmJGNSMX/e0eFpIjI/k5OeuTXrcCrhQAABaNAgDAIj4CMGgqkdLMzlaajSwxSgpEDCNX+05DFCL5xhdyZ3CTFUeioiknYqtYfCQjIznTWygrqjyd62wlEvERAKAWxEcAAIv4CMCgqUTyH6dylNQ16ilG28ltfC/hbJP3VBl5hsgWM8CF2djQ2XJcbPGcnN3NiZLkvqVtGDobAFAL4iMAgEWjAACwuKcAoAHKUz13AmKlmqGKeigHondzqvsLzvwIcoIDeR5yvTi+HHAvNh2nnGtTlqe68zp4ylMZEA8A0B3ERwAAi/gIwOApT622p7NvXgZN9naWsY8vSvJNWiDnTcgkR1JOrFQmPnJisHLTisbX06MZAFAL4iMAgEV8BGBwqjZKKtfbuSCqkmRUI6MhOReDnNZTzo8pY5xs8jmp2Dn5K4sy6Xt706MZAFAL4iMAgEV8BKBhB8pzYpsyczBIZedjKG3jRExyEDs5HaenCirW0c45lrOdf1rR6PUYEA8A0A3ERwCA3m8U1q5day5l0jzOOussZ98JEyZU3Ketra233gqA/hYllR7ebYrRo9yhiqF9ePfRz338CMPoYTqjlR66+ijhETqPon34tjePfL7mR/jxw+1m10/uKQwZMkRNmzbN+7z+pf7CCy+Y5cMPPzxxm/32209tu+22ic9lZHkWAKAmvdYojB07Vi1evNj7/O23367OPPNMNXToUPWVr3wlcZubbrpJHXXUUT14lgDQ2PpN9dHChQvNvyeddJIaOXJkX58OgMHKV/XTZbNQbOaJnWRnNxnRiAoiX4WS0/FNVhiJfm9bj+WpXopvl/DaSa9VSb/IXPT9ht///vdmWV8tAAAa+EpBR0f6Rs0uu+yijjnmGO92N998s7r++utVa2uriaOmT5+uTjvtNDVixIhePV8AGKz6vFHQjcEvf/lLs3z66aeXvWF8zz33OF/fdddd6vLLLzf/zpw5s+JrzZ8/Xy1YsCDVea1cuTLVdgAGcKe2eGVR4Omw5pvFzRdFyVnYPL/SvLFS/HegGEfJiZl8Cl2PVUV61PeNgo6NXn311bLRka5Guuyyy8yVgb6a6OjoMDetv//976tly5ap2bNnqz/84Q9qypQpZV9r3bp1aunSpT3yPgBgMMj1lxvM+hf+xIkTE7fRVwLSsGHD1AknnKCOPfZYdcQRR5hf9BdffLF69NFHy77WuHHjKjYc8kpBx1QA0EiC0PS66BubN2829wY2bdqkbr31VnX22WdXfYyHH35YzZo1y8RO77//vho1alRdzm3q1KmmsRmhRqlDghl1OSaAXtSNMZG6KDdDm93EN/tZinoe377xobPTSIjgn239D7WxuN78UVzqD+bdXfWhRYsWmQZB/+V/yimn1HSMUke3YrGoVq9eXeczBIDGkukP0dGcOXNqriBqbm62y3nTnRsAMODuKaxZs0Y9/fTT3e6bsGLFCru800471eXcADRoTNSNDm/xYbijzSsPce0bprtcpzgvOetbdJD+f6VQ6pugB7vrztAV1113nfl30qRJavz48XU8QwBoPJm+7ptwxhlnlL2Rojur6TGP9E1kSX89d+5cc19C+8EPftDDZw0Ag1+ur/om6PhINwa6USjnr3/9q7rxxhvVvHnzzFXFDjvsYEpFdcmovoegq46uueYac18CwCCLbfqj0BMHVR0reTrBlVOosqopeoH+3SiUbjB/9rOfVbvttlvZbU899VRzZbFkyRL1+uuvq+XLl6tsNqt23313deSRR6rzzz9fHXjggb105gAwuPVZo1BqGCo59NBDzQMA0AA9mgH0E40aB/WzWKmcmiIns60aWENnAwD6BxoFAIBFfAQ0gt6MhqqqiqlRNdFJXwuLdfvMaomcqsWVAgDAolEAAFjER8Bg0ctDRfep7pxff4yewp4+p/SxUz//zgMAehONAgDAolEAAFjcUwAa7d5BN/J475ST/ZC3fDPt+w/74b2HXsCVAgDAolEAAFjER8BgjYxSxCRVx0HdLVVN83p16rWbZhrMsj2Eg0xDxkpcKQAALBoFAIBFfAQM5JgoZZyTKibyHSvFvuXmWa9atvKx9GyMFRXTxGfF+lUvhb0QK6V5393ElQIAwKJRAABYxEfAIImMykZEVUZD3jgok+m7uRtEdOJ9hWIxMYbyxk1lIqbAEy2lipXqGSX1QmQkcaUAALBoFAAAFvER0BfSRCzVdj4rt32amMgXDaXZxnfMOkpVcSRf2xc3iYgpLPd2ipnaY6WeipJ6AVcKAACLRgEAYBEfAQMgMvJWFsnty1QfeWOiFNGQu2+KuKqHhtcOfFU/Mp4R5+qNm8Q2QXybYopoScZK4rXl98gbJfnOu4ZIrKdwpQAAsGgUAAAW8REwECKjFFFNl6qfFDFRqmjIWd+9+ChNZZK/o1mK+EhWHMnt00ZMgSda8sVKKSqUer2zWzdxpQAAsGgUAAAW8REwWCKjePVQNuvZzvMa2UzF+Mg5TpqKpho+my4VQSWeOEfmOU4c5IuPCqJiKB7bBJWrl4IqK5RSRUlp9UIlElcKAACLRgEAYBEfAf01MkoTE5WJjwJfHORZ724fJMZQ7utVOW5STCjfn7eyyFeJJGIfsa8TPRUK8sWixUBGT7H4qJAcLYV1qlDyRkn9qBKJKwUAgEWjAACwiI+AARwZBTLaiXcac57LJsdEzjaeDmtifeiLmHz7pvzMnHgmTec1p/pIxEci/nFiMhkxZYvJEVM8WipGzwW+CiW5v+wUWK8oaeuT7tc9XInElQIAwKJRAABYxEdAvfVmZOQsl+m8Jo8l98llK8dEzr5iG2d9cqe2sFzE5sQqyZs41T1OZZEc40jEMDI+EuudWElGT7HzC4IoDgqLno56skJJ7NsjUdLWJ3s1SuJKAQBg0SgAACziIyCtGiakTzOWUb0ioyAX+985TUyUS7NN5chIdkRzK5Rq+AydyEjsKmMbEbGEnigpyMtIRsQueVFVFD+frPg8xXZOJKbyiafdE1FSX3Ry40oBAGDRKAAALBoFAIDFPQWgHvcRfHMf+GR64D5Cl3sK4l5AU67ivQbn3oFczqZY77m/EL+nEKbpBC4rKmWcLnsly/sIheT7CEHOs17eL4kP5Feo3IvZuXfQjfsL8d7Uae5LpZqPoZvlqVwpAAAsGgUAgEV8BPR1z+V6RUaypDRebirio7BJREbiNcImGTfJWCk6v6KIj4o5T2QkqyZjg+N5p68U5Jh0zjwLnvWZfHTQjIyJ8iJi6iwklnwGsie22S6IztXT09nZvhtRkvwZ8JWnlu/F3DPlqVwpAAAsGgUAgEV8BPRkz+U0g9355iKoMjJyKozizzUnx0dFGRM1y2hILDd5IiPxckWn+sgfH8ksRVYipao4EglJRlQfFUVq40RJndEOGXHecj6FjIiLylYZBcnrnV2dt5BPjnacD8czqF/8T/W0A+clntTHZ1VFERJXCgAAi0YBAGARHwFpVdtBLc4zp4F3eswUlUjeAe3ikZFYLjaL+EhERgWxLCuO5PqijIycqiRffKS6WX0UJMdHorLIXQ4SP7KiXO9USrlhUMYTHTrRkIx9nMH7xAt65oFwIyaPlJ3aeqoSiSsFAIBFowAA6F589Pbbb6vHHntM/fGPf1TPP/+8WrZsmWptbVVTp041X5fT2dmpbrjhBnXnnXeqV155RbW0tKgDDzxQfetb31InnXRS2X1fffVVddVVV6lHHnlEvfvuu2rMmDHquOOOU9/73vfUbrvtVstbAXqv4ijwrZf7eqbQ9FUZNSVHROUio0KLJz5qEXFQk4yPxLJYL6MkXyWS6mZ8pDwVR4GIg5zIqFPMsyBitoxYn3a60Izv/MSy08HNeQ+ymkjmWOI9iJIrp7opdk41TeHZ243C3XffrS644IKq92tra1MzZ85UixcvVtlsVu27775q8+bN6sknnzSPiy++WF177bWJ+z777LOmAdi0aZMaPXq02n///dXq1avVL37xC7Vo0SLTSB188MG1vB0AQHfio5EjR6oZM2aoSy65RN17773qhz/8Yar99C993SDov+r//Oc/q+XLl5urhQceeMBcMfzoRz9SDz30UJf9tmzZok4++WTTIJx99tnqrbfeMlck69atU2eddZbauHGjeV5frQAAevlKQf9i1o+ShQsXVtznnXfeUTfffLNZvvXWW9Xee+9tn5s9e7a66KKL1D/+4z+qK6+8Un3hC19w9l2wYIFpAPbYYw/1s5/9TDU3N5v1Q4YMMcd8+umnzVXDLbfcor797W/X8paAnq84StNJzVeJ5KlccsYoEsvlIqNCi6gmakmuMso768Uxm+SypyNbrp7VRyoxJgpkhzURDTkVUSLGymWcHnHR+aT8/mZ81UQythFDeDvVZWL70CnNqtypraZqpG5WIvXajeYHH3xQdXR0mF/sRx99dJfn586da/5dunSp+QUv6XhIO/PMM22DUKK/1lcL2q9//esefAcAMPj1WqPw3HPPmX+nT5+e+PxOO+1kbxaXttUKhYK9ee3bt7Re3/jW2wMA+nnntVWrVpl/9ZWCz8SJE9WaNWvUSy+9ZNetXbvWXGGU21fvp7W3t6vXXntN7b777onbzZ8/30RRaaxcuTLVdmgcTrVRXHcqjuSyd1Y0z/ps8jhG5mvZ6Uwui2go35JcfVRoUcnrfVGSpxIpzIb+mdc8ZTxyHKSg4ImMnPhIVB9t/VXx8RO+iqPoPedkeZN57YxYlsN2Z5OroHIiGhIVR26nNk/EJM7DV4lUdojtQjdmZ+sPjcL69evNv9ttt513m9JzGzZs6LJfuX3lerlvnL4voeMpAEAfNwq6HFWL3xOQdAWSJquISvuV27e0X3zfuHHjxqkpU6akvlKgmglAo+m1RkFXCmmlKCiJjn+0oUOHdtmvtK/8Or5ffN+km9mlG9qV6I54XFUMcjV0WHP3z9S/4si3rxMlyRnSxLKIiOLVRG4ElLw+L/7XKgzxREkyPnKiJNlprEx8lEkRH8mEpZC8LCOjjPiV4ry2jFScsblFJNWlPKqYeLKyyigsiM9cdiYTEZMTE2WSv9epKpG2biiO62RGyfsPlBvNusNZPA6KKz1X2ja+7NtXrpfbAwD6aaOw1157mX91ZzWfUilqaVttwoQJNjby7VvaT8dIu+66a13PGwAaSa/FR4ceeqi67bbbTI/mJG+++aapPCpta08wlzNRjh7m4plnnlGf/exnu+yr12sHHXSQGT4D6O0Oa/Hxaioe11lOHspZjt0jIwg3PhLVR2L91u2SxzJyOqa1eCIjGSU5lUhh5fjIiZJilTBO9ZGs0JFTsnmqj0Rk5FY4iWM641Elv7BbGRSr9Ckmj6kkP9tAdhgUMVEgZn2T3zv5Hpzvte/nIVYR5a1ESjOdmj1u0P+uFE488UTV1NSkXn75ZTPOUVK5qDZ58uQupadz5syxPaf1gHqSvs+gGxvtlFNO6cF3AACDX681Cp/85CftTd5zzjnH6Yugxzu67rrrzPIVV1zRZV+939ixY018dN5559mKJP2v/lrHRzvuuKP6+te/3ltvBwAGpZriozfeeMP8RR+v/tED3G2//fZ2vR7PSD9K9C/+F154wURBeoTU/fbbzwxyV7oncOGFF5orirjhw4ebgfc+97nPmVFRf/Ob35gOanoobd0vYZtttlH33XefGjZsWC1vB6iNrzObr8Oab8Yv3xhHvsjImXlNRkRlqo/k8Nci9vF3WBPLQ8Lk+KglOTJSTSJSycU6hzkdysQTMtEREU6Y93wGzkxvlSMjt7pJrI919MrIuEpUGTmVT/mi53shthfn7Ru3SonObt7Oj+WG95bvu9DHVwp6KIn333/fPvQvdi2fzzvr9eimki4Xfeqpp8xoqJMmTTK9nN977z115JFHml/6119/vfc1p02bZhqdM844wxznT3/6k/lXj4ek18v7EACAXrxS0BVBcp7SauhKovgVRFp6OIs0I7ICAPp59REwEHlnV6ulw5p3HCRPxZGMJjxRkqyKkcNGb/1aJc+YJsYscqIkJ1byREZDRWTULGKilii/yIrIKJtzc42MiI+c8X5EFVBRxDuFvBjyW3bUaxdVV/LzluTfreKYToc4MYbS1tcLKg/JLWMiGWl1ejqpyUokOfy3bwa4eCzpG1bb05GtuzOyMUczAMCiUQAAWMRHaCzdHe+o6terrsOaEyU5VSvJEZOMNWSHrvhzMkqSHdmc4a994xq1JEdGmSFR9pJrjqKM5uZofVPWjY+yKeKjgoh6OgtRTNTREb2JvNi3KH6NFZ2xiJIjo2I++bMwX3cmR0ZhLnl9xon7kr933gq0VB3Ztr7D3sSVAgDAolEAAFjER0CNUo13VHW1ki+OkPFR8no5UX25aEkuu7FSWHEsI1llJCOjlpZo+JkhTVE+0xyLj2SclJERUOiJjOTYQs7w15FOsW8oOp8V8+K8ZWRUJnILnc8m+bPNpPheBPJ7keZ7XUPvs6rHQUqJKwUAgEWjAACwaBQAABb3FIAa51Do1iB4aQY58/R0dgaAy5SZxVF8LTNxeR/BNy+Bb4A72VtZlp7K+wjDmsT9hZw71H1O9Lb13VPIi/falhc1s/K9yR7QYuC6vDjXUPQeDjuS7yPIz8J8nZX3CMLkz9b5zD290cv1UK60TWx7b8/nHsKVAgDAolEAAFjER0Bf80QCTlTgTGPpiS9isYM/8lCJ02X6lgPPAHeyvFSWnsrIaEjWHXGuWXzti486Csm/lvKit7KMmDrFORV880447015PzPl+5x8n7Pv+yLjxL7sUV8DrhQAABaNAgDAIj4CeolvDoXYRhXXi6TFjYjKDcPvW05xLDmFppwPQQ5uJ6MkWWEk4yLzdUb2aI62K8ZLpz6WF+vla2QzueQ5GmTFUFDley7zOakUx+rO99Ttndy3uFIAAFg0CgAAi/gISMvXCambulTA2Ner8ji1nJ6nekbJweec1Ct5PoRMiuWtX4tKJuc52aktU/FYvvNI8x7iMVuPfLYZz3HkoHnVn0aZ1wu6M7aee6hunwwAYNCgUQAAWMRHQFrFWF1IbP6CWgXiuGE3ZmH0TDdQXuhbluU2YlHOXSDHH0qxvPVr+XdocvVRmmP5ziPNe6ilvCeodp9i5e91XdXxuFwpAAAsGgUAgEV8BPSSMBQVM2I5tlHF9U5RjYwpYrvK57zLKY4VFkWEI5YLxeQpNOW4RL5xjOJVRr6xj+R4R/I1Cp5zkucq30Oq91zmc1IpjpXme+dbL382+hpXCgAAi0YBAGARHwF9zRMdOBGTp3pGVrPEK1sCGaX44qOC6ExVSF4f5qO/HQv5bGKc01EoVpwtLT6WkXfmNREZyWN1iNdzoiRxTvJc3ffmWY5XCRV9UZLnc/Z9X6qNkvoRrhQAABaNAgDAIj4CtFBmBTX8reSpJHGiBid28MQInm0CEc+4x0xXSZMpiDGExIT2gRjZ2omPOuVE9yI+ahKVRR255PGHBBkFxYe/9sVHbiwVLbd15hJfuyAiI9WZSXwPznsT7znjjuyt5OfkrUTyxUrie+R8f30dy3zbxLZ3KpN6IX7iSgEAYNEoAAAs4iOgRk5MVO3OxWLlZSfKkDFF8noZfZiv82I7T2Qk45OMExlF68OcWN8uOql5IiM5FpHsyBafMU1GTnIfX6c4GRnlO7KJ5yQjo4xvWUZmsfgocD6b5M82zfdCfu9Sfa9r0FMd3rhSAABYNAoAAIv4CI1FXnL7JlSv6+sVPcuZyhVHIl4IC8lVLhlReVOMV9LI+EPER5nOMDlW6RCnkU2eLayYEa8nfn10yiGuC2K8opw7DVgmkzxjmjMMtxxfSXRMk1VGTmTUITrEtSe/H2e5Uy7HIzflid8862XFkRMlpahEcqqKPD8nfYArBQCARaMAALCIj4AyQllpkolf1ouMRcYFMpaSy74Obp5OaqHTYU0s52V8lBwlbf06eu2siEmKopoo64mMnIngPSlbUXRMC8V4Qvmm6JwKouObOZSIj5zjeobqDlN0TJORUda3LN5n1hMlmec6ZeTmi+KSvxdOZCS/d55Obd5OafHObmkqlkTkJH9ma8GVAgDAolEAAFjER0Ct5GV6NkVM5BvfpiAqdGRnL09kFOTkellJ5EYLWVlBJI5blOcq4pzQGxn5ZjMTFUPiPEIRW8nXNV/L48rOb6HnNeQw177OaDIacuIj5VkOE6OkrV+HYrmY+NnKz9wZkyqfpkOi+F57fh66dEpLM05WHXGlAACwaBQAABaNAgDA4p4CUIe5FVINjldl72anPNV3f6FTzE8gSk3NoTy9mHNOWajI/50B7pIzflWsPEdB2CHvKcTmBpAfp6ckNc10od57Cs79BbHcFiauz4n7C1uPG4r35ClDFZ+57z6CU4Zax17MVQ+CZ4+bfj+uFAAAFo0CAMAiPkLj6u7geM4lf5W9m2U0JPKSIOspZcyLyEIMSid7WTv7mupWz2uL1wudeExENU4pbXKEIwfgk72knbkYZC/pLiWp8gmx2leS6p0HwhcfhSlKUmOfWXvyc4FcdspTk79HvpJU2Utd/gy45amxKClFL+Z64koBAGDRKAAALOIjIKX4QGOBjGdkj1Txp1bgHQTPU3kie7yKyCeQVUIyPhK9k2VFztYVsrJILkf755SIQuSJy5xHVkTJyKgpOTIqyvhIfkbx6iMPp/pILPum0XTnh0iOkmSVkdNTOV591CGec3o0i/hOfo9krCfjI7kst/d8372D48Wf8/Ri7u4geBJXCgAAi0YBAGARHwF17sjm6E4lkiz1kQPLiThCDhIXr6DKpKiociMjOTeAiKVEBVBBdIgritjGiY+cAfdUt+IjkW6pjJzu0lOJJOdDkIPbyYhJVhjJiGjrcwXxnPyck5fdAe6SY8CqK4581UZpdbMqiSsFAIBFowAA6F589Pbbb6vHHntM/fGPf1TPP/+8WrZsmWptbVVTp041XyfZtGmTeuihh9Tvfvc7tWTJErV27VpVLBbVTjvtpI466ih1wQUXqP322y9x36eeekodffTRZc/pK1/5irr77rtreTtATR3ZvFN1ymkqq61EKorYQXYak1Uu8hzEclDDX3/ynOTJyqk8g4JY75nW0+m85szRoOpWfeRWIslxiZKn0HQ6n+U91Uex+CgjI6OOfOKy6hTrxbKsOJJjVcnvaZqKoy7jGzljJ9Vv2s26Ngr6l6/+JV6Nb37zm+rOO+80y0OHDlV77rmnaRRefvll9Ytf/ELdcccdav78+eqss87yHqOlpUV95jOfSXxun332qfJdAADq0iiMHDlSzZgxw/yC1o9Vq1apSy+9tOJ+xx9/vDr//PPNvs3NzWbdhg0b1Le+9S111113qXPPPVcddNBB3iuGsWPHqsWLF9dyygCAnmoUzj77bPMoWbhwYcV9brjhBvWJT3yiy/rRo0eb/ZcvX67+/Oc/q1tuucVsC/Q78aqOaquRPFUl3gjIV63kdGoTkUUs9vLFSRmnM1Q2cbrQYk50kCv4xjgS24uKI6egKdZ5zRmRW76l0DeMthyDKU0lUvL0pM7Q155OaV2qjDrylTupOZVFnqokbyWSp3KpFnUcB6nXbjQnNQglTU1N6thjjzXLL730Um+dEgCgv/ZTaGtrM/8OGzbMu81HH32k5s6dq1avXm3ip4kTJ6rZs2ermTNn9uKZAsDg1S8aBV259MADD5jl6dOne7fT9x8WLFjgrPvpT39qrjL0ze/tt9++7OvoG9nx/X1WrlyZajsMQr1RiVRtpzYlKl485xCknQ1OrJcxTCiqjIKmbGL0ImMlGSVlZExUJj6qvvpIRl3J693qo+SKI7fzmZwhzR8fKVll5BvXKB9tE4rlajuplR3fqDsVR9XO1NZfGoXvfe976p133lE77LCDc6+iRFcrnX766eq0005TkyZNUmPGjFHr1q1T9957r7riiivU448/bq4YnnnmGZXNxgZwF/Q+S5cu7eF3AwADV583Cvov/J/85Cdm+ec//7mpbIo75JBDzEPadddd1YUXXqgOP/xwc3Xx7LPPmmPphsNn3LhxasqUKamvFPQVDAA0kj5tFB599FF1xhlnmOWrr75anXjiiVUf47DDDlNz5sxR99xzj7r//vvLNgr6foR+pKE74nFVgVrGRepWpzY5rpHvZVNESV2OK9fL8xMVRKE8PxldiXGXAhEfZcT60Kk+knFY7NxTpHGpKpFk9VEheSY073rfbGldxpXKi/XF6iIjp+KoULmTmqeD2tYvw8aYee3pp59WX/ziF1VHR4e65JJLUvVz8NFXC5ruLwEAGGCNgo56dEe2LVu2qO985zvqmmuu6dbxSh3h8rLVBgD0//johRdeULNmzTJjIekezPXoqLZixQrz7/jx4+twhkANlUjVDrEtIwGnV1amflFSKDqmyfUyGhIVMIHsyJYX5yEio7BTvDdPZCTjpi4nlaaayxkTSCXHQfLzk1FS0RMTyaofX0e0MtspGQd1IzJyq49qGLsoTWRUQ8VRn10p/O///q867rjj1Icffqi+9rWvqZtvvtmZZrAWb731lh1TSR8bADAAGgU98J3uZLZ+/Xr15S9/Wd12220qI/4qKkdv/8gjj3SJh/Roq3ocJd3I7LjjjqlvIgMA6hgfvfHGG2ry5Mn26/b2dvOvHr9IdiC76KKLzEPTg97pvgjaa6+9po488khv2eiiRYucdbpB0Ot0fwXdi3n48OHmCkGfRyk20sNyJ5WzAnUTvyxPcZXrrUSS1Uqyk1mxPlFS/Hyd44q+PM7Q2b6KIxElBbLiSJxfIGdbE59LfDwmh/ybsJhmaO/KkZHbaSy5EskbEZn9C9HLeWZS605k1GVY7OjFxGItsVLYt41CoVBQ77//fpf1+i95uV7fSI43HNr//M//eI+t+x/EXXvttaZa6cUXXzSNgR7uYptttjFVR7rTmr5CGDVqVC1vBQDQ3UZhwoQJ/hbPQ0+UU6vzzjvPPAAAg7xHMzDo1KtTW7VRUrwyRXaUyngiI7ks46CCHMtIdlIT28hoSN4fFJVIXcKjaquPJN8YQjIm88VEvn1jYx+Fzv6eobDFsaqOjGqZRa2HOqn5MEczAMCiUQAAWMRHQE92bOvFKCk+LrVzRjLykNvJDmFOp7MguVopLyMjT3yUMi6SHd6czmipoiRPNORUKHnGHJIRUTyaKSQ/58ZKxZ6NjNLGRXWsOJK4UgAAWDQKAACL+AgYCFGSHARIVtvI14rFCc5w2XI8at95OHlVJrkTmG+WOBkflaswEs95t/J28KpcfZRqOOpy8VHRcyyn+ijljGkDKDKSuFIAAFg0CgAAi0YBAGBxTwEYCPcXnIxf9lQWm8icPf5yzr2HTFX3GmTpqXMcZ96IQnKpakyaofK9Q+ikyePT3FPw3TfobrnpAL6PIHGlAACwaBQAABbxETDQoiTJ0wPakAPneUpXfbGSnODAHfjOFx+pVPFRt4KQNDFMuWioUkQUey5ME2N5IqBuRUa9HBfFcaUAALBoFAAAFvERMMCiJMnXA7rLYTPVxUrOvikiJuec0syZUINUE3uVi4YStil7zGKdps4cAJGRxJUCAMCiUQAAWMRHwECIkqRqK5TKREupYqUqYw4nbqqnauOj7nSI66nKon4aGUlcKQAALBoFAIBFfAQMhCipGxVKWzfzRBtiyCLnWLKbmeyMJucVkLv2VGTUrUqk7sU8Yb0qiwZAZCRxpQAAsGgUAAAW8RHQn/jihWorlFJ3fpMd06qLmJzTqGVUI/naaaKaNKqMc1JFRIM0JvLhSgEAYNEoAAAs4iNgMFYo1dD5LY2yHeSqlVzI1C3Vvp+EA6ja9x24kZHElQIAwKJRAABYxEfAYK9QqiUi8VUZ1atKqLd1JxYapDGRD1cKAACLRgEAYBEfAYNFuVij2rGJuhO3eKKnfhUBpXq9UDUirhQAABaNAgDAIj4CGkGaKKRew1/3dsxTrQaNhdLiSgEAYNEoAAAs4iMAtcUqvTzbmhdxUF1xpQAAsGgUAAAW8RGA2hDbDEpcKQAALBoFAIBFowAAsGgUAAAWjQIAgEYBANAVVwoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAIDuNQpvv/22uvPOO9W8efPUtGnT1LBhw1QQBOozn/lM2f0mTJhgtiv3aGtr8+7/6quvqrPPPluNHz9etbS0qJ133lmdc845as2aNbW8DQBAPSbZufvuu9UFF1ygarXffvupbbfdNvG5TCa5nXr22WfVcccdpzZt2qRGjx6t9t9/f7V69Wr1i1/8Qi1atEg99thj6uCDD675nAAANTYKI0eOVDNmzDBXBvqxatUqdemll6be/6abblJHHXVU6u23bNmiTj75ZNMg6CuFf/mXf1FDhgwxVxXnn3++uu2228zz+jyGDh3K9xUAejM+0r+YH330UXXNNdeYX8bjxo1TPWnBggVq3bp1ao899lA/+9nPTIOg6X9vvvlmNXHiRPXXv/5V3XLLLT16HgAw2A2IG806HtLOPPNM1dzc7Dynvz7rrLPM8q9//es+OT8AaOj4qLv0X/fXX3+9am1tVWPHjlXTp09Xp512mhoxYkSXbQuFgnr++efNst4uSWn9H//4R7N9Npvt4XcAAINTnzQK99xzj/P1XXfdpS6//HLz78yZM53n1q5dqzo6Osyyjo+S6PhIa29vV6+99prafffdE7ebP3++iaLSWLlyZartAGAw6dVG4fDDD1eXXXaZ+ct+l112Mb/sFy9erL7//e+rZcuWqdmzZ6s//OEPasqUKXaf9evX2+Xtttsu8bhy/YYNG7yvr+9LLF26tG7vBwAGm15tFPSVgKT7N5xwwgnq2GOPVUcccYT5hX3xxRebm9glst9C/H5Cie6zUKIjKR99Q1w2OJWuFModCwAGoz6Jj+J0GenVV1+tZs2apZ544gn1wQcfqFGjRpnnSpVGmr6ykF+X6NhIHstn7ty55pHG1KlTuaoA0HD6TfWRjpa0YrFoOqWV6I5qSVGSJNfL7QEAA7RRkNFQPp93hsYoPffKK68k7ltqRHSMtOuuu/b4uQLAYNVvGoUVK1bY5Z122sku53I5E+VozzzzTOK+pfUHHXQQ5agAMBgaheuuu878O2nSJDPgnTRnzhzz78KFC1VnZ6fznL7PoIe50E455ZReO18AGIx6rVHQndX0mEfvv/++s15/rW/+lnot/+AHP+iyr35ed3LT8dF5551nK5L0v/prHR/tuOOO6utf/3ovvRsAGJxqqj5644031OTJk7tU/yxfvlxtv/32dv1FF11kHpoem+jGG280w23r+wQ77LCDKfnUpZ/6HoIeHVWPpVS6KpCGDx+u7r33XvW5z33OjIr6m9/8xnRQ00Np634J22yzjbrvvvtMiSsAoJcbBT2URPwvfk3/cpfr9eimJaeeeqoKw1AtWbJEvf7666YB0cNR6F/uRx55pBnt9MADD/S+pp63Qe/zj//4j6Yfw5/+9CfTsJx44ommN7SvFzMAoIcbBf2Xvv4FX41DDz3UPLpDD2eh7ysAAAb5jWYAQN+jUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAACwaBQAABaNAgDAolEAAFg0CgAAi0YBAGDRKAAALBoFAACNAgCgK64UAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAACwaBQAABaNAgDAolEAAFg0CgAAGgUAQFdcKQAALBoFAIBFowAAsGgUAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAACwaBQAABaNAgDAolEAAFg0CgCA7jUKb7/9trrzzjvVvHnz1LRp09SwYcNUEATqM5/5jHefK6+80myT5vH73//e2fepp56quM+pp55ay1sBAAg5VYO7775bXXDBBVXts8suu5gGxOf1119Xb7zxhho6dKiaPHly4jYtLS3ehmefffap6nwAAHVqFEaOHKlmzJhhfkHrx6pVq9Sll15adp+zzz7bPHyOPvpo0yicdNJJ5vhJxo4dqxYvXlzLKQMAeqpRiP+CX7hwoeqOtWvX2sjozDPP7NaxAAAD/Ebz7bffrsIwVDvvvLM65phj+vp0AKBh1XSlUE+6MfjlL39pls844wyVyfjbqY8++kjNnTtXrV69WjU3N6uJEyeq2bNnq5kzZ/biGQPA4NXnjcLTTz+tXn31VdsolLNhwwa1YMECZ91Pf/pTdeyxx5qb39tvv33Z/efPn99lf5+VK1em2g4ABpM+bxRuu+028+8RRxyh9thjj8RtdEXS6aefrk477TQ1adIkNWbMGLVu3Tp17733qiuuuEI9/vjj5orhmWeeUdls1vtaep+lS5f22HsBgIGuTxuFzZs3q/vuu88sn3XWWd7tDjnkEPOQdt11V3XhhReqww8/XE2fPl09++yz5mpBNxw+48aNU1OmTEl9pdDa2pr6vQDAYNCnjcKiRYvUpk2bTOe3U045paZjHHbYYWrOnDnqnnvuUffff3/ZRkHfj9CPNKZOncpVBYCG06fVR6VS1pNPPlmNGDGi5uPoqwVN95cAAAzARkH3TdA3mStFR2noSiQtn8/X5dwAoFFl+vIqQZejTpgwQR111FHdOtaKFSvMv+PHj6/T2QFAY8r0dd8EXVWkB7Sr1VtvvWUG59OOO+64up0jADSiPmkU9JAWa9asMY1Bpb4J2pe//GX1yCOPdImHlixZYsZg+vDDD9WOO+6Y+iYyAKCO1Ud64Do5kml7e7v5d/ny5U4Hsosuusg8fDeYP/vZz6rdd9+94uvpBkFXKun+CroX8/Dhw80Vgj6PUmz00EMPeQfSAwD0YKNQKBTU+++/32W9/ktert+yZUvZvglpB7+79tprzU3pF1980TQGeriLbbbZxlQd6U5r+gph1KhRtbwVAEB3GwV9c1jfF6iF/it/48aNVe1z3nnnmQcAoAFGSQUA9A80CgAAi0YBAGDRKAAALBoFAIBFowAAsGgUAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAACwaBQAABaNAgDAolEAAFg0CgAAi0YBAGDRKAAALBoFAIBFowAAsGgUAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAYNEoAAAsGgUAgEWjAACwaBQAABaNAgDAolEAAFg0CgAAi0YBAGDRKAAALBoFAIBFowAAsGgUAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMAKwjAMoy9Rst1226kNGzaojMqq4WoEHwyAAWuz2qiKqqBGjx6t1q9fX3bbXK+d1QDT1tZm/tUf5Eb1QV+fDgDU7fdaOTQKHmPGjFF/+9vfVEdHhyoUCmro0KHqU5/6VPe/Kw1g5cqVqrW1lc+Mz4yfs35izZo1pkHQv9cqoVHwWLt2rfl36tSpaunSpaZBeOGFF+r7nRqk+Mz4zPg5G7i40QwAsGgUAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFsNcVPCNb3xDrVu3To0bN67SpuAzqxk/Z3xm/QVDZwMALOIjAIBFowAAsGgUAAAWjQIAwKJR8HjyySfVCSecoHbYYQczg9g+++yjLr/8crV582bVaPQ03v/93/+tLrnkEnXEEUeoT3ziE6qpqcl8Nscdd5z6t3/7N7NNkgkTJqggCMo+0kwROBBdeeWVFd/7zTffnLhvZ2en+vGPf6wOOOAANXz4cDNn+DHHHKPuv/9+NZgntqr0eZUeZ511lrNvI/+c1RslqQluuukmNW/ePPOLbvz48WrnnXdWf/nLX9RVV12l7rvvPrV48WLzP2mjeOKJJ9SMGTPs17vvvrvabbfdzBR/jz76qHn86le/Mp9NS0tL4jH2228/te222yY+l8kM7r9N9BSIe+65Z+JzSaXO+pfXzJkzzc9ZNptV++67r/ljRP+hoh8XX3yxuvbaa9VgM2TIEDVt2jTv8/pzKc1+ePjhhydu08g/Z3UTwvH888+HmUwmDIIgnD9/flgsFs36N998M5w6dar+czg86aSTGupTe/TRR8PddtstvPHGG8N33nnHee6Xv/xl2NLSYj6Xiy++uMu+u+66q3nuySefDBvNFVdcYd77GWecUdV+3/nOd8x++jP/v//7P7v+gQcesJ/1gw8+GDaahQsXmvc+dOjQ8MMPP3Sea+Sfs3qjUYg58cQTzQ/X6aef3uXDWrVqlWkw9PPLly8PG4X+H7Cjo8P7/NVXX20+k+222y4sFArOc438P2stjcLbb78dNjc3m/2eeOKJLs9ffvnl5rkpU6aEjeaoo44y7/20007r8lwj/5zVG9dTwqZNm9TDDz9se5jG6QhA57raokWLVKMYOXKkuYfgM2vWLPPv+vXr1bvvvtuLZzb4PPjgg6qjo0Ptscce6uijj+7y/Ny5c82/S5cuVatXr1aNQt9v+P3vf2+WzzzzzL4+nUGNewrCsmXLVHt7u8nFDz744MQPbPr06eqxxx5Tzz33XG99j/o9eQNP35RPom+oXn/99aq1tVWNHTvWfI6nnXaaGjFihBrsli9frr761a+qt99+27zfT3/60+rUU0819wriSj9X+vNJstNOO9n7OXrbiRMnqkZw++23m3t8u+yyi/3DLEkj/5zVTd2vPQawW265xVyC7rnnnt5t7rzzTrPNzjvv3Kvn1p99+9vfNp/JAQcc4L2sT3psv/324SOPPBIO9vgo6aHvWX33u98N8/m8s88RRxxhnteRnM+MGTPMNjpKagT6vt7uu+9u3vP3vve9xG0a+ees3oiPBB1/aOUqi0rPbdiwoX4t8wCmY4xSWaUuWY3TVSILFixQK1euNBU0+nN76KGH1OTJk9V7772nZs+ebY4xGOm/VC+66CLzF72O1fQV1Z/+9Cd13nnnmb96b7jhBnXppZc6+/Az2JWOjV599dWy0VEj/5zVXd2bmQHs//2//2f+spg+fbp3m8cff9xsk81mw0anb4russsu5vP40pe+VNW+W7ZsMTdL9b76L99G86Mf/ci891wuF65Zs8auL/1FfOutt3r3/drXvma2Oeecc8JGoG/UV/r/0qfRf85qwZVCrE5a0zf6fPQ9h3LZeaP48MMPzQ3m119/XU2dOlUtXLiwqv3153f11VfbfhAffPCBaiQXXnih2nHHHVU+nzd/0ZbwM+jSf/Xr/i+13mBu9J+zWtAoCKNHj3Yu4ZOUnitt26hVWp///OfNjXl9s/R3v/udqVCqVqkDUrFYbKhKGk13SjvkkEPM8qpVq+x6fgZduspP/7wNGzZMnXLKKTV91o38c1YLGgVhr732Mv/qv371MANJSj9UpW0bzZYtW9Txxx9vcnL9GehKLD3sRS2am5vtsv6LudGU3r9876Wfq1deecW7XyP9DJauQOfMmVNzBVGj/5xVi0ZBmDJlivkB0hHRkiVLEj+wZ555xvx72GGHqUajb5SeeOKJ6umnnzZjzTz++OPmZmqtVqxY4ZRaNprS+9dDqZQceuih5l89xEWSN99805Sjym0HK/0+9c9ad/smNPrPWdVquhMxiH3hC19I1aP5xRdfDBuJ7tH8d3/3d+a9jx8/Pnz11Ve7fcxTTjnFHG/SpElho/mP//gPWzK5ZMkS5+Z9U1NTxR7NkydPDge7UknvhAkT7HAztWjkn7Na0CjE6P9BdQ15fOyjt956y4599MUvfjFsJLqWvvQ/1tixY03jmMaPf/zj8J//+Z/D9957z1mvv/7GN75hfykuWrQoHGxWrFhh3mP8jwc9DMhdd90Vjhw50rz3448/vsu+3/rWtxLHPtLjHZXGPvr3f//3cDDT/9/p96/fq24cymnkn7OeQKOQ4J/+6Z9Mo1DqpKb/Kiv9z7j33nuH7777bthI9C+x0v9Y+q+2adOmeR9Lly61+82bN8921NL/gx988MHh/vvvb8ow9Xp91aVLMwejZcuW2c9Mjwmlf4YOOuigcPTo0Xa9LrHcsGFDYhnlYYcdZkufdafAiRMn2v0uvPDCcLDTYxiVfnYqXZU28s9ZT6BR8HjsscfCWbNmmf+hdYOw1157hZdeemm4cePGsNHcdttt3t6i8YcckOzZZ581I34eeuih4Y477mg+x2HDhpnP8txzzzW/OAcr/cv+qquuMlcCuu/BiBEjTCz0yU9+0vxc3XHHHV16M0vt7e3mF5n+5aZHBd12223DI488Mrz33nvDRuqboN9zJY38c9YTAv2f6u9EAAAGI6qPAAAWjQIAwKJRAABYNAoAAItGAQBg0SgAACwaBQCARaMAALBoFAAAFo0CAMCiUQAAWDQKAACLRgEAoEr+P7uxzW7IQNyjAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data = tobac.testing.make_sample_data_2D_3blobs(data_type=\"xarray\")\n", "\n", "plt.figure(figsize=(8, 8))\n", "plt.imshow(data[50])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 14, "id": "ba44cd13", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:30.918586Z", "iopub.status.busy": "2026-02-02T20:12:30.918475Z", "iopub.status.idle": "2026-02-02T20:12:31.072141Z", "shell.execute_reply": "2026-02-02T20:12:31.071492Z" } }, "outputs": [], "source": [ "%%capture\n", "threshold = 9\n", "dxy, dt = tobac.utils.get_spacings(data)\n", "features = tobac.feature_detection_multithreshold(data, dxy, threshold)" ] }, { "cell_type": "code", "execution_count": 15, "id": "91ed6612", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.074103Z", "iopub.status.busy": "2026-02-02T20:12:31.073934Z", "iopub.status.idle": "2026-02-02T20:12:31.078909Z", "shell.execute_reply": "2026-02-02T20:12:31.078266Z" } }, "outputs": [ { "data": { "text/plain": [ "60 501.0\n", "61 30.0\n", "62 325.0\n", "Name: num, dtype: float64" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features.where(features[\"frame\"] == 50)[\"num\"].dropna()" ] }, { "cell_type": "markdown", "id": "68e3b85d", "metadata": {}, "source": [ "Obviously, the feature with only 30 datapoints is the rightmost feature that has almost left the imaging area. If we now use an *n_min-threshold* above or equal to 30, we will not detect this small feature:" ] }, { "cell_type": "code", "execution_count": 16, "id": "95284244", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.080804Z", "iopub.status.busy": "2026-02-02T20:12:31.080682Z", "iopub.status.idle": "2026-02-02T20:12:31.192353Z", "shell.execute_reply": "2026-02-02T20:12:31.191791Z" } }, "outputs": [], "source": [ "%%capture\n", "features = tobac.feature_detection_multithreshold(\n", " data, dxy, threshold, n_min_threshold=30\n", ")" ] }, { "cell_type": "code", "execution_count": 17, "id": "caeea658", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.193914Z", "iopub.status.busy": "2026-02-02T20:12:31.193815Z", "iopub.status.idle": "2026-02-02T20:12:31.197820Z", "shell.execute_reply": "2026-02-02T20:12:31.197451Z" } }, "outputs": [ { "data": { "text/plain": [ "60 501.0\n", "61 325.0\n", "Name: num, dtype: float64" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "features.where(features[\"frame\"] == 50)[\"num\"].dropna()" ] }, { "cell_type": "markdown", "id": "72563272", "metadata": {}, "source": [ "Noisy data can be cleaned using this method by ignoring smaller features of no significance or, as here, features that leave the detection range." ] }, { "cell_type": "markdown", "id": "31c70494", "metadata": {}, "source": [ "## Minimum Object Pair Distance `min_distance`" ] }, { "cell_type": "markdown", "id": "7bf8c298", "metadata": {}, "source": [ "Another way of getting rid of this specific feature is the *min_distance* parameter. It sets a minimal distance of our features. Lets detect the features as before:" ] }, { "cell_type": "code", "execution_count": 18, "id": "b2a7de25", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.199116Z", "iopub.status.busy": "2026-02-02T20:12:31.199034Z", "iopub.status.idle": "2026-02-02T20:12:31.330719Z", "shell.execute_reply": "2026-02-02T20:12:31.330250Z" } }, "outputs": [], "source": [ "%%capture\n", "\n", "data = tobac.testing.make_sample_data_2D_3blobs(data_type=\"xarray\")\n", "\n", "threshold = 9\n", "dxy, dt = tobac.utils.get_spacings(data)\n", "features = tobac.feature_detection_multithreshold(data, dxy, threshold)" ] }, { "cell_type": "markdown", "id": "e6fd198e", "metadata": {}, "source": [ "A quick look at frame 50 tells us the found features and their indices:" ] }, { "cell_type": "code", "execution_count": 19, "id": "f7db8149", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.332313Z", "iopub.status.busy": "2026-02-02T20:12:31.332218Z", "iopub.status.idle": "2026-02-02T20:12:31.336372Z", "shell.execute_reply": "2026-02-02T20:12:31.336021Z" } }, "outputs": [ { "data": { "text/plain": [ "array([60, 61, 62])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 50\n", "mask = features[\"frame\"] == n\n", "features_frame = features.where(mask).dropna()\n", "features_frame.index.values" ] }, { "cell_type": "markdown", "id": "6a24c009", "metadata": {}, "source": [ "Notice that *to_dataframe()* was used to convert the Dataset to a pandas dataframe, which is required to use the *calculate_distance* function of the analysis module. The distances bewteen our features are:" ] }, { "cell_type": "code", "execution_count": 20, "id": "b15da9d2", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.337787Z", "iopub.status.busy": "2026-02-02T20:12:31.337701Z", "iopub.status.idle": "2026-02-02T20:12:31.340465Z", "shell.execute_reply": "2026-02-02T20:12:31.340061Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(77307.67873610597)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tobac.analysis.calculate_distance(\n", " features_frame.loc[60], features_frame.loc[61], method_distance=\"xy\"\n", ")" ] }, { "cell_type": "code", "execution_count": 21, "id": "4afc4eb3", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.341600Z", "iopub.status.busy": "2026-02-02T20:12:31.341525Z", "iopub.status.idle": "2026-02-02T20:12:31.344066Z", "shell.execute_reply": "2026-02-02T20:12:31.343646Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(64289.48419281164)" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tobac.analysis.calculate_distance(\n", " features_frame.loc[62], features_frame.loc[61], method_distance=\"xy\"\n", ")" ] }, { "cell_type": "code", "execution_count": 22, "id": "9a5ff68a", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.345223Z", "iopub.status.busy": "2026-02-02T20:12:31.345145Z", "iopub.status.idle": "2026-02-02T20:12:31.347783Z", "shell.execute_reply": "2026-02-02T20:12:31.347307Z" } }, "outputs": [ { "data": { "text/plain": [ "np.float64(101607.57596570993)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tobac.analysis.calculate_distance(\n", " features_frame.loc[62], features_frame.loc[60], method_distance=\"xy\"\n", ")" ] }, { "cell_type": "markdown", "id": "3ee0e593", "metadata": {}, "source": [ "With this knowledge we can set reasonable values for *min_distance* to exclude the small feature:" ] }, { "cell_type": "code", "execution_count": 23, "id": "50dd4ad8", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.348857Z", "iopub.status.busy": "2026-02-02T20:12:31.348779Z", "iopub.status.idle": "2026-02-02T20:12:31.350321Z", "shell.execute_reply": "2026-02-02T20:12:31.349998Z" } }, "outputs": [], "source": [ "min_distance = 70000" ] }, { "cell_type": "markdown", "id": "15c4336f", "metadata": {}, "source": [ "and perform the feature detection as usual:" ] }, { "cell_type": "code", "execution_count": 24, "id": "ffc101e8", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.351360Z", "iopub.status.busy": "2026-02-02T20:12:31.351284Z", "iopub.status.idle": "2026-02-02T20:12:31.472046Z", "shell.execute_reply": "2026-02-02T20:12:31.471510Z" } }, "outputs": [], "source": [ "%%capture\n", "\n", "data = tobac.testing.make_sample_data_2D_3blobs(data_type=\"xarray\")\n", "\n", "thresholds = [10]\n", "dxy, dt = tobac.utils.get_spacings(data)\n", "features = tobac.feature_detection_multithreshold(\n", " data, dxy, thresholds, min_distance=min_distance\n", ")\n", "\n", "n = 50\n", "mask = features[\"frame\"] == n" ] }, { "cell_type": "markdown", "id": "3e0d7a87", "metadata": {}, "source": [ "Plotting the result shows us that we now exclude the expected feature." ] }, { "cell_type": "code", "execution_count": 25, "id": "1d1d0131", "metadata": { "execution": { "iopub.execute_input": "2026-02-02T20:12:31.473360Z", "iopub.status.busy": "2026-02-02T20:12:31.473277Z", "iopub.status.idle": "2026-02-02T20:12:31.524977Z", "shell.execute_reply": "2026-02-02T20:12:31.524401Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAKrCAYAAAAebZGTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAASARJREFUeJzt3QmYFEWa//G3ae5LUUEuORUdxAPaA0UEFXBYFRwFx5VHBR3F8dFB1+cvPt6zI+P5OLrOroAXiuuFOqvu7I6CqICLiwLi4rCiHJ7gBR7IDfV/Ip6p2OjsiuzIzKqs6qrv53nKzq7KzMoqm3orfhkRWZXJZDICAICINOJdAABkURQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAA23KLz22mty6qmnSvv27aVFixZy0EEHyQ033CA//fRTsQ8NABq8qoY099F9990nkyZNEnXIXbt21YXhr3/9q2zbtk1+9rOfyYIFC2SvvfYq9mECQIPVYFoKixcvliuuuEIvT5s2TT755BNZsmSJrF69WmpqamTFihVy0UUXFfswAaBBazAthdNPP11eeOEFOe+88+TRRx+t9diHH36oY6Tdu3fLsmXL5NBDDy3acQJAQ9YgWgqbNm2Sv/zlL3r54osvrvP4AQccICeeeKJenjVrVurHBwDlorE0AEuXLtXnDZo1ayZHHXVUznUGDx4sc+bMkbfeeisvz9mjRw/56quvpHnz5tKzZ8+87BMAimHNmjWydetW6dChg6xdu7bhF4WVK1fqn926dZMmTZrkXKd379765wcffODcjzoXMX36dK/nVOcsVLK2ZcsW2bhxY6zjBoBSor7o1qdBFIUNGzbon2E9i7KPhX2Ar1u3Tp+cjkJ1e1U9m9KQ5PSOz7a++7fXy9d+o762qqqqSOuHbZOv++Oul9Z+ABfVEUd9wVXJR1kUBdXsUZo2bepcR0VLinrhLp06dZIBAwZEehNVQVA9n6J+qIWt73rM9WEcdVmdcM+17Fon+PuuXbtyLrv26/McSYuC/bu93KjR/50Wq66uznl/1HVcy8Fjsh9z3e86bp/XlkbBRGWoqanRX4h9ovAGURSy1W379u3OddQ5h+w3e5eJEyfqW5Q3Me43bJ/78/nh71MUXB/2wd+jFgjX/WFFKElRiPrhby/bx+Fax/V+u4pA8Pjs53BtY6/v+vvggxzF0CB6H7Vr165WjJRL9rHsugCAMi0Kffr0MSd/d+zYkXOdVatW1VoXABBdg4iP1HkAdT5BRUSLFi2SQYMG1Vln/vz5+ucxxxyT6rHFOTlcCpHRzp07ax2T/btrPZ99RT3v4Jutu2Kixo0b1xsHudax77ePz77fl+s8gs/fhysmsrd17dM3YoqzDSpTg2gptG7dWk4++WS9nKtLqRrRPHfuXL08ZsyY1I8PAMpFgygKipoJVX3DmTlzpi4M2W8+qpvp3//93+tvemoqjMMOO6zYhwoADVaDiI+UI488Uu6++275h3/4B92D6JZbbpF99tnHzJJ64IEHygMPPCClwrevf76iJFcPIPscjL1O8NyM6zE7GnLd7xM3+bwfvvGRKwJyLdsDHu37fWI5ezlprBQVMQ+KocG0FBQ1S+rs2bNl5MiR+voJqiB0795drr32WnnnnXd0kQAAVEBLIeukk07SNwBA/jW4olAK4gxSi7q9T88i17Ir5nFFQcHf7UGC9rJPrORaxyeecY1UDouDXNGQfb8dY9n3J42GXNsUoveRz/7piYSKi48AAIVFUQAAGMRHBRLW28anJ07UHkeumMhedsVCYY9l55QK3p+dpDC4jitKcvWO8hmgFoxq7IkR7eXspIiKPRukz/xN9ra+XL2lfF6fjbmPUEpoKQAADIoCAMAgPkqpx5FrmyQ9jnziI59YKPi7HQ3Z16ewl+11XFGS65h84pJg1GL3GnJFRva06fZrbdmyZaQYK+k1Hlzr+Mxl5LOt7zo+vZGYEwlBtBQAAAZFAQBgUBQAAAbnFPLI51xB2DY+10rwuaaBq1uo67xB8Hc1r1R9y65zCq5zGD6T44V1SXWdU7C7ntqvz+c8gitP97kedNhj9rL9fEnOETA5HtJCSwEAYFAUAAAG8VERFGsUc1h8tGnTpnqX7fho8+bN9e436ohmOyKxr5kQjI/sbqj289ldUl3XdfB5bp9YKOwxn0tz5utSm2HrEz8hDloKAACDogAAMIiPPIRNaFeoy3FGHdFsL/tMaGfHP8GY6IcffshLTySfayvEiY/sHkf2sus5ok7AZz+3vewbH9nbuP6f+vRESnItBiAuWgoAAIOiAAAwiI8SinpthLBtovY+8umJZEc49oR2dvwT/N1edvVEsuMne9kVH9kxln3cNjt2CV7q0u5xFHWAnM2OeezniLoc/N31/8h+TVF7IqWtVI4DxUVLAQBgUBQAAAbxUREU4noKrvmOXJfQtKOk4O92HOQTK7liKft+n5jHjnbs+Y2CA9Ncg9R85hOy4xyfazTY69jLwdfhiozs++3XF7UHW9SBbL7bAEG0FAAABkUBAGAQHxWI72U6o8ZHrsjCFSUljY9cUZJ9vz3Y7ccff8y5T/s47MjH5op2gsdu78sVsbh6MrnmULKP1b7fXg4et2tOJdf03FGX8zmNNj2L4IuWAgDAoCgAAAziI09Re4uE3Z9k6mz7ftfgNZ95kOzlYDxjP2bHKvayHRlFjY9cg9fs3jnB+Mg1DbctamRkv2Z7DiX7WMOe1/We+/y/8+ETJdnobYR8oKUAADAoCgAAg/goBt9psZPsK8mgNtfcR67l4O+uuYV8BrjZ97viI9dx2/FPcPCaa/prOyayoyF7e1ds5vPeuCKi4DH5REY+y0Cx0VIAABgUBQCAQXxUIL5XW/NZpxBRUjAKcUVGrumvXQPhXL2Yog5ec/UwCq7n6k3kEwf5xESu5aSRkQuDzFBsFAUURYdMRibs3i3H794trdUEeyIyf/duebRxY/mKidyAoqEoIFXNMxm5Y+dOOXf3bql9Kllk2M6dct3OnTKzulquqa6WbRQHIHUUhSIrRJTkMydPMMLxmV/JFbf49nZSBeH5XbtkSMhrVoXiwl275IDNm+X0pk1l698Kg8/U4K51XK8n6pxSwf9Xrh5RUQcn+iBWQlo40YzU3P23glDfx6F6XMVKdwa6zQIoPIoCUrFvJiPn/e3bcX1zfGYfP3fXLn3uAUB6iI9KiE+8kK+4KWw+pqhTePsM4hq/a1edcwj1Ueufu2OH3FFd7dW7x+eYCjWYLMkANAavoZTQUkAqjo/5oTnEkdsDKAyKAlLRJuZ2qrsqgPQQH5Uo1zTIPtMj2+u4luNsb09t7bNf+/4f1XKM1sKmqir9vD7Hbh+fvew6Jtf9caagZtpqlAtaCkjFvJhjDuY5PtwBFAb/4pCKR6qq5P9GFPhR6z9qTWkBoPAoCkjFl1VV8tjfWgs+4xSUx6qrmfICSBnnFIosaZZd3358zg8Ef/dZtiels5ftS1/ao4rVdQ/+X3W19Nm+vd6eSOoo5zdqJJObNZPqvx2z61KbrmXXsfq8Nt9zMlHP+yQ518M5C6SFlgJSo6asGN20qTwQEiWp+x9s1EhOb96cuY+AIqClgNQLw6XV1XKzmiU1k9HTXrTJZHTvJHVSeUajRjoyasJkeEBRUBSKIGoU4Io2fKIdV7xixzxhj7kud+laDrt8pX1834vIPSLyBytKyr6Gpjkux2k/h318rktw2suu98AVgflEUsHfo8Z3QKkiPgIAGBQFAIBBfBSDKwZwzXnvO+FZWO+gXPcniY/s5WDcYscz9v0tWrQwyy1btsx5TYOwyCjXOvZ74+rFFIyJWrdunfM47GVXlORadkVmYfFR1NHeUUeBA8WYVJGWAgDAoCgAAAziozxyRUZhg56iDjSz73f1nrGX7YikefPmOeOY4Hr2Y3ZkZA9Ga9OmTc44yD5WO3rZtm1bvdc6sLcNxkd2NNSqVaucUZJ9rPZrsF+3Kz5y9WJy9VYK/h51skAbURIyJXQxKVoKAACDogAASD8+Us2jhQsXyosvvigLFiyQFStWyA8//CB77rmn9O/fX84//3w555xzcjale/ToIR9//HHo/rds2VIrJsi3OL2Jou7XJz7yiYxcvYdckUrwd7s3kT0YzV52vQeuQXB29GTvx2a/zuDgNTsasuMje9nVE8ne1l52xWn2cliPrSRzLRXiehm+26Byo6GSKgpz586VYcOGmd979eolPXv2lDVr1sjs2bP17cknn5TnnnuuTt6d1a9fP9ljjz1yPubqwgkAKNGWgioCV1xxhZx99tnSoUMH89jMmTPloosukj//+c9y0003yW233ZZzH/fdd58MHTo0rUMGgIqTWlE46qij5IMPPqjTq0Q599xz5dNPP5XrrrtOHnjgAfn973/f4L/5h8VNPtFB1EFqrqjGp1dR2AA0u6eQvWxz9Rqyn8O1/ziD11wRkCtKcsVHrnVc8ywF4yOf3kdxpuTOhem1S1+mgcVELql98rZt2zZnQcgaOXKk/rlhwwb5+uuv0zosAEApjlPYunWr81ts1tSpU+Wuu+7SJ5U7duwogwcPlnHjxtXqMw8AKIOioE4yK4cddphuVeTy9NNP1/r9iSeekBtuuEH/HD58eL3PMW3aNJk+fbrX8ajeUfkSbPq7BrZFnRbbNceR3SKze/rYxTbYA8g15XXUyMg1WM7ev2uf9usPRjWu6bJd0ZCrh5IrYrKP1X4u17JvTzDmQSo/mRKJiQp1HCVRFJYsWaJbAco111xT5/Fjjz1Wn29QLYNu3brpfFp1a73xxhtl6dKlMmrUKHnzzTdlwIABoc+zbt06/VwAgBItCl9++aX84he/0CdH1U/VMylItQRs6hveqaeeKieddJIcd9xx+oN+8uTJultrmE6dOtVbOOyWgoqpAKCSFLUofP/99/oE8yeffCI1NTUyY8aMSNuruGDKlCl6H2ocxHfffacHw7lMnDhR33yo48m2KlQE4Gqq+fQKCet95LPsmnbaFR/ZEY69flgs5OoR5DOIzhXtuAbBJY2PXNN5R42P7GV7W58ptcMiI58Ba0mm1C7UQDQGuNWWZkxUKpGUUrR+n5s2bZKf//znOv45+OCD5eWXX3aeSwijoqXsB82qVasKcKQAUDmKUhQ2b94sp5xyirz11lvSp08fmTNnjuy9996x9mV/q3NNnwAAKNH4SHU9HT16tMybN0/PafTqq6/q7qVxLV++3Cx36dJF0hB1HqSkvY9cA7zsGMaONuz7XZFR8Lh9Ih3XtNh2ZGRPke2a78jnuYLTVLsiHbvXUJLeR665j1zTaAd/TxIfuSRdhzioNGKcTIGioULtN9WWgvqQOPPMM3XLoGvXrvo8gPqZxB133KF/9u3bN/G+AKDSpVYU1DdWNdDsP/7jP3TLQBUENRdSfdRgNTXn0bffflvrfvW7Omk8a9Ys/ftvf/vbgh07AFSK1OKjZ555xnyAqyb/hAkTnOuqIqCm01Y+++wzuffee2XSpEk6bmrfvr3uKqq6jKpIQjXPb731VhkzZowUm6u5HmfuI5+eSK44yL7fjj/sdcKanq6YyNXLyB6Nbvc4suMjV+8mn6gqeByuAXI+cyJFnTo7Tu8j1+C1fM2DhNKMYDIJ9lXo44iy/9SKgp01r127Vt/CuqpmqXEL6gUtWrRId11dtmyZ/kenpt4eMmSIXHrppXL44YcX/PgBoBKkVhTGjx+vb1ENHDhQ3wAAFTCiuaFw9Tjy6YkUFglE7Ynk2tZetmMO32ZjkkFqdivQFR/5XLUtrPeRaw4iV+zjGozmWt8191HYldeSREYuacdKlRJXlUoPoIzH+sWOtBr2RQsAAHlFUQAAGMRHBZL0guo+kZEdWSRtckadnts1x5EdGdnLUa/gFtb7yDUPkitistdxLbu2DYuP7N8LMd+Ri0/vNRQ3hsnkMSYqxHOHoaUAADAoCgAAg/go5Z5IYfvyWccVK0XtcRR8Xp8rvdnPYfcmsqMkn6m6fXofBY/PJ9JyRUw+0ZB9v6vHVbBHlOs98+lxVIhYyVc5x0xpR0aZBOtEvT/uelHRUgAAGBQFAIBBUQAAGJxTSEGcDNfnkpgudg7u2+XV55yCvV9Xd1PXdRN8JuMLy9Ndx+SavM9n2efcgWs5+NxRzxHk69oK5Xx+oBgyKZ4jiHNOIY3LdtJSAAAYFAUAgEF8VOTuqT7xjs8IYBc7LgmLoXwuu2l3MbXjFlfXU9c1FAoVH/l0pfWJhnz2H/ZYvmKifEZG5RwzpXEdg4xH1JMkJsrnccR53EZLAQBgUBQAAAbxUQp8r6fg2san6e8TJYXt0xUf2XGQK57xuSxo1Ca+7/H5xEc+k/3FGZ2cZu8jpBMZZQrQgyhJxJSP7aOipQAAMCgKAACD+MiDat5HnWQu6YR4rp5MUfdj32/HOb69e+weRK7rN7iulZCkqesbHyWJlaLuM+yY8hUZ0eOouL2MCtHjKGmUlI8eTvQ+AgDEQnwEADCIjzxFjYbC1vcZ8BZ1v/ayHXn4Njdd0YarJ5IrPrIVqvdRktgnyTphx1gqkRGSzWOUKXAcFGc/+YyifNBSAAAYFAUAgEF8lIJgDOAzR1K+eiXZwmIln2NybZ/PuV1c8tXTp1CDzEohMirHuKkQU0XHmWcoEzG28Zk23negZ5xtgvfT+wgAEAvxEQDAID6KIc4gNZ/tk/RKSvK8wd/TmM8lqiQxTJJoJ2l85HOsSZRjZJRE0rmLMgl6CkWNeXx78/kMDq3v3yPxEQAgFuIjAIBBfNSAoySfiCnOFeAKcZWpfMpXj54kUVKhjslHuUdGUf+GkqwfpyeSzScy8omJ7DnG7OU4+yU+AgDkDfERAMAgPmpgUVIaA+ribJ9mlJSvqcRd98eJfJIck49yj4wKIel02ZmIA9Nc0Y5r2RUZBXsfuR6LMhcZvY8AALEQHwEADOKjBhYl+ayTNGpIMh+Taz9J1vdViB5AREbl3eMo6eC13QliojjxkWu9+gbFER8BAGIhPgIAGMRHKShUT59CT7WdNMby2WccUbcvVHyUz23ysS2S93jLFGC+I5/IaOfOnTnv993eXs51TK6rI+ZCSwEAYFAUAAAG8VED7JkUZT++PYAKMT132PPF3U8cpRgTFWI/yN+V05L2PooaGYXFRzt27Ii031zHR+8jAEAsxEcAAIOiAAAwOKdQhucX8jl6OF8jpZOc/0iqFLuIch4hf6OYk/y7CRvRbPO5jGaS8wj2sn0OIewx1zkMRjQDAPKG+AgAYBAflZB8defM5z6jxhylOKK5WPss5H6RW5JLyfpuk/EYuRy1q6odC/nGR1Em12NEMwAgFuIjAIBBfFRBPZRc+wxTzJHI+ZL2MZXie1CJ8nk5zoxHfBS1J5IdBbl6Iinbt2+vdz1XrER8BABIhPgIAGAQHzUwheihFOf5kkh78FqhNLTjRelejnOHZ+8jO0qyl+sb1MaEeACAWIiPAADFKQo333yzbnKH3aZOnZpzW9U8uvPOO+Wwww6TVq1ayV577SUnnniiPP/882m+hJLlej9LUTkca6keb7lRsUf25nN/1P3EOY7du3ebm+t+103FObluYeuo+Ch7U5FRnFuU112UcwodOnSQAw44IOdjnTp1qnPf1q1bZfjw4bJgwQKprq6Wgw8+WH766Sd57bXX9G3y5Mly2223pXDkAFDeilIURo4cKTNmzPBeX33oq4LQs2dP+c///E858MAD9f0vvviinHXWWXL77bfLoEGD5LTTTivgUQNA+Sv5cwpffvmliZQeeughUxCUUaNGydVXX22iKdRVX1xXyVEI70v6ksY4hTyeNG67HPFR2E1F51FudmyUva+seh+p1oB6cfvvv7+ccMIJdR6fOHGi/rlkyRJZtWpVEY4QAMpHUeKjZcuWyTnnnCPr16+XNm3ayKGHHipnn322PlcQ9NZbb+mfgwcPzrmvLl266FhpzZo1et3evXsX/PgBoFwVpSi8++67+ma3BqZMmSKTJk2Su+66S59Mzlq5cqX+qVoKLqoQqKLwwQcfhD7vtGnTZPr06V7HuGLFCqkkhZgiO22VGoNVukJckc1X1EFtUafgDv7uM0DOHuyWfY6S7X3UsWNHfQ7gjDPO0B/kqpWgPvT/5V/+RZ83uOeee6Rp06b6xHHWhg0b9E/VBdUl+9jGjRtDn3/dunU6ZgIAlEBRuOSSS+rcd8ghh8j999+vIyDVy+juu++WX//619KjRw/THVVRxcKlWbNm+ueWLVtCn191dx0wYIB3S6G+/QFAuSmZuY+uuuoquffee+WLL76Ql156SS6//HJ9f/PmzevM8xG0bds2/bNFixahz6FOSmdPTNenpqaGVkUIohqUonxGRhmPuY+S7McVKwWvkuYzj1J9U3g3yN5H6jzC0UcfXes8gtKuXbtaMVIu2cey6wIA4imZomBHRPaJkj59+uifH330kXO7bFfU7LoAgDIoCsuXL9c/u3btau4bOHCg/qlGNOfy+eef655H9roAkIaMY5BavvaTr1uDLAp//vOf5f3339fLI0aMMPePHj1amjRpIh9++KGe5yhXN1Olf//+od1WAQAlVBTUB746yasGrtnUCZInn3xSD2ZTTjnlFDnyyCPN4/vuu685OXzhhRfWGougTkjfcccdevmmm25K6ZUAQPlKrfeRmn9DDRxTNzWuoHv37tK4cWN9riA7vkCNWn788cfrbKs++BcvXiwLFy7Uo5779esnmzZtMucSVM8l1aIAgGL1wqtKMHgybD+ux6Iul1xLQY07uOWWW3RLYM8999TFQI1qVieX1aypM2fO1PGQeixIdTV9/fXX9aC2vn376t5J33zzjQwZMkSeffZZPQoaANCAWgrqw/66666Lvb0qHmo0dHZWVABAGQ9eA4Ck7Lgk6UC2qjxFMq5tGzVqVO/9wd/teeHsZddAtuw6UY65ZHofAQCKj6IAADCIjwCUZWSUzygpSTTkiolc69uxUPB3136D2wSjJOIjAEAsxEcAAIP4CEBq0oh0oijUQLEqx7Kr95BrWVHT/OTqZeR6/3L1SiI+AgDEQnwEADAoCgAAg3MKAMq+u2nScxlVEUci+3QddXU1DZ5TUBOHZrmuVW8fR67Ld3JOAQAQC/ERAMAgPgJQUcKilKqIXUyjRka5JqsLdjsNcsVdrufINTke8REAIBbiIwCAQXwEoCzF6XFUlSA+co1ctp/bjo/syMheJ+xYfSbRyxUfBa/REIaWAgDAoCgAAAziIwAVNdjN5hvVVEW8DoK9X591wnofuY7d53Kc2ecgPgIAxEJ8BAAwiI8AFEW+5i9KMq9R2KCujOM57CjGXscVDSU9Jtdz25GR6zoLDF4DACRCfAQAMIiPAJSlKPP9ROmxVOWIc2yuKClpfOSKjFzxUXaZuY8AALEQHwEADOIjAGUjn1dYK8SgNh/BSMpnSm7X/EoMXgMAJEJ8BAAwiI8AlM1ANp/1k6ry6H3ks61r2Y5/wp4vV0yUa/vgPupDSwEAYFAUAAAG8RGAspT0yms211TYdlTjEyX5xEf2QLSw53PFR7leK/ERACAW4iMAgEF8BKCkFLMnUsaxfdQrurmiJJ/9B7eNGhkRHwEA8ob4CABgEB8BKBtRr85WqEFtux29klzPba8fPCbX9q4oycbU2QCARIiPAAAG8RGAsp0KO8k+qyJGS64Bbj69hMKOybWN6zlybcvgNQBALMRHAACDogAAMDinAKAs+Zw78B3dXBXxMp1Rz3/EOafgWsfn/jC0FAAABkUBAGAQHwFoEAoVz0TdpipPo6B9958kPnLtMwwtBQCAQVEAABjERwAqSliUVJUgMvJZP58xluu5k6KlAAAwKAoAAIP4CEDFTpQX1tOnyhEBpdkTKc72cR630VIAABgUBQBA+kVh7dq1ugnjc5swYUKtbXv06FHvNlu3bk3rpQAoIfbnQBr7qnKs4/v5luSmrqGQvUXZriTPKTRv3lwGDRrkfFx9qC9evFgvH3vssTnX6devn+yxxx45H7MvOAEAiCe1otCxY0dZsGCB8/FHH31Uxo8fLy1atJBf/vKXOde57777ZOjQoQU8SgCobCXT+2jGjBn65xlnnCFt27Yt9uEAqHBVCabeTtJDKaxHVNTnjqMkMhd1vuGNN97Qy6q1AACo4JaCio5UdevWrZuceOKJzvWmTp0qd911l2zZskXHUYMHD5Zx48ZJmzZtUj1eAChXRS8Kqhg89thjevm8884LPWH89NNP1/r9iSeekBtuuEH/HD58eL3PNW3aNJk+fbrXca1YscJrPQDlN6gt6X7zFT2FrZfkuUu6KKjYaPXq1aHRkeqNdN111+mWgWpNbN++XZ+0vvHGG2Xp0qUyatQoefPNN2XAgAGhz7Vu3TpZsmRJQV4HAJSDxqVygll94Pfu3TvnOqolYGvZsqWceuqpctJJJ8lxxx2nP+gnT54ss2fPDn2uTp061Vs47JaCiqkAoJIUtSj89NNP8txzz8U+way6r06ZMkVGjhwpc+fOle+++0723HNP5/oTJ07UNx81NTW0KgAk7vWTr2m3fSXdV1F7H82aNUs2bdqkv/mPHTs21j6yA912794tq1atyvMRAkBlaVQK0dGYMWNi9yBq2rSpWd65c2fejg0AKlHR4qM1a9bIvHnzEo9NWL58uVnu0qVLXo4NQMORz15GSSQZTOYb88TdV4OYOjs7NkFNdpdk6oo77rhD/+zbt6907do1j0cIAJWnUbHHJpx//vmhVUwNVlNzHn377be17le/q5PG6ryE8tvf/rbARw0A5a9xscYmqPhIFQNVFMJ89tlncu+998qkSZN0q6J9+/a6q6jqMqrOIajBbrfeeqs+LwGg8mKbUleVxzmK8nWlt5IrCtkTzMcff7z07NkzdN2zzz5bv3mLFi2STz75RJYtWybV1dXSq1cvGTJkiFx66aVy+OGHp3TkAFDeilYUsoWhPgMHDtQ3AEAFjGgGUBqIg9IVJwpK4/9RSUydDQAoDRQFAIBBfARUAKKh8lCVQu8jWgoAAIOiAAAwiI+AMkFEhHygpQAAMCgKAACDogAAMDinADQwxTx30JDOW6TRfbMc0VIAABgUBQCAQXwENACFim3SjoPyeWnKUnquckJLAQBgUBQAAAbxEVBCSjEmKsWIyYdPNBT2XFUlGC0xIR4AIFXERwAAg/gIKLJ8xSVx9hN1m1IfvGbHK0l7H2Uc26cdK6X9fLQUAAAGRQEAYBAfARUQGfmsl691Sp1vxFTliG1KJVYqFFoKAACDogAAMIiPgJQUegBZ0iip1COmNOKZKo9oySdWyuexFmq/LrQUAAAGRQEAYBAfAWUaGUWNifJ1f6H4xDmFilqqIj5H2pFPPtFSAAAYFAUAgEF8BJRJZBS8P0kclM/jyJeoMUzaUVIakVYasRQtBQCAQVEAABjER0CJRkZJ708SE+VrnUJFIfb6rm3T7hlU5REZNYReSbQUAAAGRQEAYBAfAWUSGYXFR1G391mO+lxJRY2MfJaDMhHjnSSRUdIoqVBRFC0FAIBBUQAAGMRHQEoKHRn5xkeFWN69e3ek1xB8LGoPokaNGsWOjMKipKoEMUy5REm0FAAABkUBAGAQHwGekg7QSrJ+nJ5BrnjHJwLK1/2+ry9JZJSv+8PWK1aUFHwsDbQUAAAGRQEAYFAUAAAG5xSAEIUanVuI0cbBLD9q/p+vZd/zHC4+3UftvD/qsn0MwfMGVR5dY6Oea0ja9TTqNkm7p9JSAAAYFAUAgEF8BORZkstUJomMgvuPGvvs2rUr57JrHZ/4KKx7qg+fbqX2cnV1db3328vB96yRIxqy77dfU5IoyVZK11mgpQAAMCgKAACD+Ago4MjlQlzHwDeq8YmGkixHjafCXlO+ehm5YiLXMdnr5Pq9kFGSzXdEc5KeSL5oKQAADIoCAMAgPgJSUujrGASjGlfss3PnzpzLPuv73B92TD69kaLGRI0bN673fp+eUnE0cgyKS4LeRwCAkkF8BABIFh+tX79e5syZI2+//ba88847snTpUtmyZYvU1NTo38Ps2LFD7rnnHnn88cflo48+kmbNmsnhhx8ul112mZxxxhmh265evVpuueUWeeWVV+Trr7+WDh06yIgRI+T666+Xnj17xnkpQGo9jqLuM+ocRXaE4xsZ+Sxv3749dvQUNqAu6oA1V88iOyayl5s2bZrz/kLNZ2Vz9UTK56U2CzXgLVZReOqpp+TKK6+MvN3WrVtl+PDhsmDBAv0/9eCDD5affvpJXnvtNX2bPHmy3HbbbTm3XbhwoS4AmzZtknbt2skhhxwiq1atkocfflhmzZqli9RRRx0V5+UAAJLER23btpVhw4bJNddcI88++6z8/ve/99pOfeirgqC+1b///vuybNky3Vp44YUXdIvh9ttvl5deeqnOdps3b5YzzzxTF4QLLrhAvvjiC90iWbdunUyYMEF+/PFH/bhqrQAAUm4pqA9mdcuaMWNGvdt8+eWXMnXqVL380EMPyYEHHmgeGzVqlFx99dXyu9/9Tm6++WY57bTTam07ffp0XQD2339/uf/++02zsHnz5nqf8+bN062GBx98UC6//PI4LwkoiDiX0Yyyflh85Ip6VIRb3/12ZOQTH7n2E9YjysU1f1GTJk3qjYxcPY7sKMnnvQ/jGmjnMzAv6fOV1YnmF198Uf9xqQ/2E044oc7jEydO1D+XLFmiP+BtKh5Sxo8fX+t/rqJ+V60F5ZlnningKwCA8pdaUXjrrbf0z8GDB+d8vEuXLuZkcXbd7DeL7Mlr17bZ+9WJb59vIgCAIg9eW7lypf6pWgouvXv3ljVr1sgHH3xg7lu7dq1pvrq2Vdsp27Ztk48//lh69eqVc71p06bpKMrHihUrvNZD5Qhr+heix1HUnki+8ZGrZ1G+ll3xVNhrckUydhxkx0d2YuBajhrXBZ/b5pqPKWqU5NpnKfVESq0obNiwQf/ca6+9nOtkH9u4cWOd7cK2te+3tw1S5yVUPAUAKHJRUN1RleA5AZvqgaTYvYiy24Vtm90uuG1Qp06dZMCAAd4tBXozAag0qRUF1VMo2MwMUvGP0qJFizrbZbe1fw9uF9w218ns7Ant+qiBeLQqyluhBjEVosdR1KulBX/36WVk/zuyl+0vZq5lV3wUdkyuyMM1SM0VH9mfCfayzxxHwailKuIU3q71faKdUu2JlNqJZjXgLBgHBWUfy64bXHZta99vrw8AKNGi0KdPH/1TDVZzyXZFza6r9OjRw3wrcG2b3U7FSN27d8/rcQNAJUktPho4cKA88sgjekRzLp9//rnueZRd1xxg48Y6ylHTXMyfP1+OP/74Otuq+5UjjzzS68pJQLHla7rsOPGRq9eQHRnZ59PUjAK5IiN7Hdc+g8fkinRcUY3979nVy8h33qVczxWMYxpFvNKbazlqbFjs6bKL0lIYPXq0zgQ//PBDPc9Rru6iSv/+/et0PR0zZowZOW3nldk/QFVslLFjxxbwFQBA+UutKOy7777mJO+FF15YayyCmu/ojjvu0Ms33XRTnW3Vdh07dtTx0SWXXGK+raif6ncVH3Xu3Fl+9atfpfVyAKAsxYqPPv30U/2NPtjsVBPc7bPPPuZ+NZ+RumWpD/7FixfrKEjNkNqvXz89yV32nMBVV12lWxRBrVq10hPvnXzyyXpW1D/96U96gJqaSluNS2jdurU899xz0rJlyzgvB4glau+RfM1xFGfwms+gMzsasiMje9nVXdyOnlxzJcWJj+zBa65BcWFXeqtv/3Hio2or0nI9t0+UlHTwWaGmzo7VUlB/fN9++625qQ/27P8o+377jynbXfT111/Xs6H27dtXj3L+5ptvZMiQIfpD/6677nI+56BBg3TROf/88/V+3nvvPf1TzYek7rfPQwAAUmwpqB5Bcft4qxNEwRaELzWdhc+MrACAEu99BDREvl9+4syzE3c/SXsfuWIl18A0OzJSF8Wqb4CbHfMEO4ZEjY/sAWuuQXGuSTBdkY8dSQV7K+70uJqcvWxvX4i/gTCuyChprMQ1mgEABkUBAGAQH6GipHHRdp/njjoPkk/vo2DPm6hTZ7siI1cvI9c6+ex9ZEdGrmmxfa7a5hoEZ8dTwePdaS27rujmWo7a06yUBrLRUgAAGBQFAIBBfAQUmU8vo6hRUliPmaiD2nym0XYNXgv2Poo6dbYd7/jELa5eRq7XZl+LJXh8u6zlqJGRz0C2UkVLAQBgUBQAAAZFAQBgcE4BKCCfrok+97uya1eOHfzdzsddI45dl+l0nSPwOQeRdESzzwR3rnMHrtHQYSOum1rdVX3OHYS9/7nWcd0fdpnOtLuo0lIAABgUBQCAQXwElKgkI52Dv7u28ZkAzl7HdU0DVyQTHNGs9rvH1q1y4urV8rOvv5YWO3bIliZN5K/t28trvXrJ982bOy+p6xqh7Ho+n9cTfM8yEbv+Ro0BGwKKAoBUNNm5Uy5YvFiGrlkjjQMfmod9+aWMff99eb1nT5l5xBGyg2utFw1FAUAqBWHyG29I36++cn8YZTIybPVq6bxpk9w+ZIjssE4gIz2860ADEHUCPd/LdtrbRB3B69pnrvtVC0EVBPVsYX1p1ONqvXMXL5YHjjiiVmSUr+NzveZCTXDX0HCiGUBB7bFli46MlPo6V2YfV+urcw9IH0UBQEGdmOMcQn3U+iesXl2wY4Ib8RHgqVDxQL56sMQ5Pp9rPCRdVr2M4uj79dfy7wU+Pl+ZPP2/SONvKOlgN1oKAApKdTtNczskQ1EAUFBqHEKa2yEZ4iPAU7BZnq8owNXcjxoDxIkN7G0KtbyifXs9DiEqNaCt0MfnqypP/y8KNY9RPvdLSwFAQc3t2VN2RvzQUuurEc5IH0UBQEF936KFHqms1Ne2yj6u1ldTXiB9xEdAA+ATfwTvtwd+ueYNcl3K0rWtzz5z3T+zpkaPVA4b0ayPR8VGHTroqS7Utj7PEfX4XK85bPsqj/c87SmuC4WWAoCCU1NWqKkr5vTq5YyS1P3q8duHDmXuoyKipQAgtcKgpq54pl8/PTCtb8RZUpEOigJQonwiC1fcEfzdtY191TL7w9hettdxXc3Mvt91BbLsc29u0kT+fNhhtQamZY+vSeD47P3ay65jivp6gu9Zlcd7GzVWamiIjwAABkUBAGAQHwEF5Loge9TYwRVZuOKO4O92fOKKZOxl+wL2rquq2eu4BvIFj8k1hbWrR5DrmOzlqOu4loPP3cgjMgp7/3Otk+T+tNBSAAAYFAUAgEF8BBSZTxwRtVeMb+8bV9xiL9sXunfFP67XE+xe6rqamU/vKPuYmlujne3lZs2aRXpt9v6Dx1sdcbCcz/+7YkdDPmgpAAAMigIAwCA+QkVx9QZK+7mj3h91rp9gNOJadkUydkxkRz7BC93nOlZ7/3ZvpTjxkaunkCsysu9v0aJFzvvD4qPGjvcp6vsfda6kUoqVaCkAAAyKAgDAID4C8hA3+TT/feIqn3jBNSgt2NPH1eMoSXzk2r+9vH379pw9l/LZ+yhqZBQnPqp2vD5Xr6R8RUPB9aPuK2kURUsBAGBQFAAABvERkFJPJp95kHwiI5+ri4VFL64Ix3W/zfXcrqmsg72VosZH9mtyDUBzRUYtW7bMeb9rP8HfGzuipKi9j+JcNc9HoXos0VIAABgUBQCAQVEAABicUwCKIMmIV1fmHuw6av/uOr9g5/pRJ7hzdUO1RzEHzylEvZ6Cz7UVXF1Po06Ul0aX1KoSHcVso6UAADAoCgAAg/gIFatQk+NFjQV8uqfaMYW9vm985IqJfCIjnxHGrvgoOKLZ9dyu1+3q6upzrQSf5eDlOBsniI+iXkPBN1ZKO2aipQAAMCgKAACD+AjwFGzGu+KPqFFU1HjBJ0oKOw6fyMi17OoNZPfuseOjYKTluh6DzRXP2M8X9VoRrsjI93oK1Qkio3z2OMrniGgXWgoAAIOiAAAwiI+AAipETySfgWxxjs9nAj67N5G9bMdCrsFxYddmsLkiGVcPIJ/4yCduCnusOk/XUGgIg9poKQAADIoCACBZfLR+/XqZM2eOvP322/LOO+/I0qVLZcuWLVJTU6N/z2XTpk3y0ksvycsvvyyLFi2StWvX6qZkly5dZOjQoXLllVdKv379cm77+uuvywknnBB6TL/85S/lqaeeivNygFi9h3yuj5Cvnkh2ZFEoPvGRHRO54iM7IgobUOcS9doRrpjHZyBanMtxNkpwPYVS7XGUuCioD1/1IR7Fr3/9a3n88cfNRTEOOOAA/Qfy4YcfysMPPywzZ86UadOmyYQJE5z7UN3ejjjiiJyPHXTQQRFfBQAgL0Whbdu2MmzYMP0BrW4rV66Ua6+9tt7tTjnlFLn00kv1ttmTQBs3bpTLLrtMnnjiCbnooovkyCOPdLYYOnbsKAsWLIhzyACAQhWFCy64QN+yZsyYUe8299xzj+y999517m/Xrp3eftmyZfL+++/Lgw8+qNcFyk2+mvuumCef+/KJj+yoxRUZBeOiqHMf+RyTK1aKuhz2WKM8DVgr1R5HRTnRnKsg2CMMTzrpJL38wQcfpHVIAIBSHaewdevWOhfbDvrhhx9k4sSJsmrVKh0/9e7dW0aNGiXDhw9P8UgBoHyVRFFQPZdeeOEFvTx48GDneur8w/Tp02vd98c//lG3MtTJ73322Sf0edSJ7OD2LitWrPBaD+UnjZ5ISeICO76w4xnfKClqPGPHKD69jOzX6dPbKIwrqvE5Vp9YKfieNYrYs6gQV14L+9uI+ncT5++sJIrC9ddfL19++aW0b9++1rmKLNVb6bzzzpNx48ZJ3759pUOHDrJu3Tp59tln5aabbpJXX31Vtxjmz58fOqpTbbNkyZICvxoAaLiKXhTUN/y7775bLz/wwAO6Z1PQ0UcfrW+27t27y1VXXSXHHnusbl0sXLhQ70sVDpdOnTrJgAEDvFsKqgUDAJWkqEVh9uzZcv755+vlKVOmyOjRoyPv45hjjpExY8bI008/Lc8//3xoUVDnI9TNhxqIR6sCcSQZ1Ba1ue+KkoKPuZ7Dp9ePvV/Xld58rvKW6/d8HF++lsOeoypPkZFLPiOjBjvNxbx58+T000/Xl/G75pprvMY5uKjWgqLGSwAAGlhRUFGPGsi2efNm+c1vfiO33nprov1lB8IFrwcLACjx+Gjx4sUycuRIPReSGsGcj4Fqy5cv1z+7du2ahyMEJC/zGIVxRQJJeigFoxCfK69FjY9cPYt8exxFjY9sUeOcqPen3bOoqkCxUNL9ptpS+J//+R8ZMWKEfP/993LuuefK1KlTE7+AL774wsyppPYNAGgARUFNfKcGmW3YsEHOOusseeSRR7z7Vav1X3nllTrxkJptVc2jpIpM586dvU8iAwDyGB99+umn0r9/f/P7tm3b9E81f5E9gOzqq6/WN0VNeqfGIigff/yxDBkyxNltdNasWbXuUwVB3afGK6hRzK1atdItBHUc2dhITcudqzsrkC/BVm3UKCTJ+okjAesLmM/8Q/b9Ptu6IqOw1+zz+qJGSUljnqoCzF+UxoCzfEZRsYqCGtX47bff1rlffZO371cnkoOFQ/nv//5v577V+IOg2267TfdWevfdd3UxUNNdtG7dWvc6UoPWVAthzz33jPNSAABJi0KPHj0in3RTF8qJ65JLLtE3AECZj2gGKkWS+ZHyGSW59ht12eaKmPIpSe+eOAPLqgoQGaV9FbU4uEYzAMCgKAAADOIjoAgD29KOkqLOu+QTGfmsEyZJJJavgWJhvY9spRgZFSpyoqUAADAoCgAAg/gIaABRUqEkmc7bJ1bKpyQxTJweR7Z8DVIr1cjIRksBAGBQFAAABkUBAGBwTgFoAOcXfNbJZ96cr9eQhnx2Ba0qQLfShnAewUZLAQBgUBQAAAbxEVCGUVLwsSSixlVpd0nN5zpVKTxHIbfNB1oKAACDogAAMIiPgJQUokePb++jNEdKFzP+yGe0U5VCz6J87ycfaCkAAAyKAgDAID4CGliUFGc/SeOnhiLJdRkKuU0upfp+01IAABgUBQCAQXwElGGUFBT1WgmubUs18mgocw5VNYD3j5YCAMCgKAAADOIjoIQUaj6hJPstVORRCrFUoZ63qgHERC60FAAABkUBAGAQHwEV1EMpbL8+ivncxdpnJURGNloKAACDogAAMIiPgAYm7Sue+Tx3Jakq8/eAlgIAwKAoAAAM4iOgTMSZ+wjx3s9yRksBAGBQFAAABvERUAGiTpFdzio1FvJFSwEAYFAUAAAG8RGAos+DlARxUH7RUgAAGBQFAIBBfAQgFmKb8kRLAQBgUBQAAAZFAQBgUBQAAAZFAQBAUQAA1EVLAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQCQrCisX79eHn/8cZk0aZIMGjRIWrZsqedWP+KII0K369Gjh14v7LZ161bn9qtXr5YLLrhAunbtKs2aNZP99ttPLrzwQlmzZk2clwEAyMdFdp566im58sorJa5+/frJHnvskfOxRo1y16mFCxfKiBEjZNOmTdKuXTs55JBDZNWqVfLwww/LrFmzZM6cOXLUUUfFPiYAQMyi0LZtWxk2bJhuGajbypUr5dprr/Xe/r777pOhQ4d6r79582Y588wzdUFQLYV//ud/lubNm+tWxaWXXiqPPPKIflwdR4sWLfj/CgBpxkfqg3n27Nly66236g/jTp06SSFNnz5d1q1bJ/vvv7/cf//9uiAo6ufUqVOld+/e8tlnn8mDDz5Y0OMAgHLXIE40q3hIGT9+vDRt2rTWY+r3CRMm6OVnnnmmKMcHABUdHyWlvt3fddddsmXLFunYsaMMHjxYxo0bJ23atKmz7q5du+Sdd97Ry2q9XLL3v/3223r96urqAr8CAChPRSkKTz/9dK3fn3jiCbnhhhv0z+HDh9d6bO3atbJ9+3a9rOKjXFR8pGzbtk0+/vhj6dWrV871pk2bpqMoHytWrPBaDwDKSapF4dhjj5XrrrtOf7Pv1q2b/rBfsGCB3HjjjbJ06VIZNWqUvPnmmzJgwACzzYYNG8zyXnvtlXO/9v0bN250Pr86L7FkyZK8vR4AKDepFgXVErCp8Q2nnnqqnHTSSXLcccfpD+zJkyfrk9hZ9riF4PmELDVmIUtFUi7qhLhdcOprKYTtCwDKUVHioyDVjXTKlCkycuRImTt3rnz33Xey55576seyPY0U1bKwf89SsZG9L5eJEyfqm4+amhpaFQAqTsn0PlLRkrJ79249KC1LDVTLFSXZ7Pvt9QEADbQo2NHQzp07a02NkX3so48+yrlttoioGKl79+4FP1YAKFclUxSWL19ulrt06WKWGzdurKMcZf78+Tm3zd5/5JFH0h0VAMqhKNxxxx36Z9++ffWEd7YxY8bonzNmzJAdO3bUekydZ1DTXChjx45N7XgBoBylVhTUYDU159G3335b6371uzr5mx21/Nvf/rbOtupxNchNxUeXXHKJ6ZGkfqrfVXzUuXNn+dWvfpXSqwGA8hSr99Gnn34q/fv3r9P7Z9myZbLPPvuY+6+++mp9U9TcRPfee6+ebludJ2jfvr3u8qm6fqpzCGp2VDWXUrZVYGvVqpU8++yzcvLJJ+tZUf/0pz/pAWpqKm01LqF169by3HPP6S6uAICUi4KaSiL4jV9RH+72/Wp206yzzz5bMpmMLFq0SD755BNdQNR0FOrDfciQIXq208MPP9z5nOq6DWqb3/3ud3ocw3vvvacLy+jRo/VoaNcoZgBAgYuC+qavPuCjGDhwoL4loaazUOcVAABlfqIZAFB8FAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAAEUBAFAXLQUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAABQFAAAddFSAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAYFAUAgEFRAAAkKwrr16+Xxx9/XCZNmiSDBg2Sli1bSlVVlRxxxBHObW6++Wa9js/tjTfeqLXt66+/Xu82Z599dpyXAgCwNJYYnnrqKbnyyisjbdOtWzddQFw++eQT+fTTT6VFixbSv3//nOs0a9bMWXgOOuigSMcDAMhTUWjbtq0MGzZMf0Cr28qVK+Xaa68N3eaCCy7QN5cTTjhBF4UzzjhD7z+Xjh07yoIFC+IcMgCgUEUh+AE/Y8YMSWLt2rUmMho/fnyifQEAGviJ5kcffVQymYzst99+cuKJJxb7cACgYsVqKeSTKgaPPfaYXj7//POlUSN3nfrhhx9k4sSJsmrVKmnatKn07t1bRo0aJcOHD0/xiAGgfBW9KMybN09Wr15tikKYjRs3yvTp02vd98c//lFOOukkffJ7n332Cd1+2rRpdbZ3WbFihdd6AFBOil4UHnnkEf3zuOOOk/333z/nOqpH0nnnnSfjxo2Tvn37SocOHWTdunXy7LPPyk033SSvvvqqbjHMnz9fqqurnc+ltlmyZEnBXgsANHRFLQo//fSTPPfcc3p5woQJzvWOPvpofbN1795drrrqKjn22GNl8ODBsnDhQt1aUIXDpVOnTjJgwADvlsKWLVu8XwsAlIOiFoVZs2bJpk2b9OC3sWPHxtrHMcccI2PGjJGnn35ann/++dCioM5HqJuPmpoaWhUAKk5Rex9lu7KeeeaZ0qZNm9j7Ua0FRY2XAAA0wKKgxiaok8z1RUc+VE8kZefOnXk5NgCoVI2K2UpQ3VF79OghQ4cOTbSv5cuX659du3bN09EBQGVqVOyxCapXkZrQLq4vvvhCT86njBgxIm/HCACVqChFQU1psWbNGl0M6huboJx11lnyyiuv1ImHFi1apOdg+v7776Vz587eJ5EBAHnsfaQmrrNnMt22bZv+uWzZsloDyK6++mp9c51gPv7446VXr171Pp8qCKqnkhqvoEYxt2rVSrcQ1HFkY6OXXnrJOZEeAKCARWHXrl3y7bff1rlffZO379+8eXPo2ATfye9uu+02fVL63Xff1cVATXfRunVr3etIDVpTLYQ999wzzksBACQtCurksDovEIf6lv/jjz9G2uaSSy7RNwBABcySCgAoDRQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBBUQAAGBQFAIBRlclkMv/3K7L22msv2bhxo7Ro0UJ+9rOf8cYAaLBWrFghW7ZskXbt2smGDRtC16UoOLRs2VK/iQBQLtSX3M2bN4eu0zi1o2lgOnToIF999ZVs375ddu3aRYshxrcSWlm8Z4XE35m/NWvWyNatW/XnWn0oCg5r167VP2tqamTJkiU6Qlq8eHGE/w2Vi/eM94y/s4aLE80AAIOiAAAwKAoAAIOiAAAwKAoAAIOiAAAwKAoAAIOiAAAwKAoAAIOiAAAwmOaiHhdffLGsW7dOOnXqVN+q4D2Ljb8z3rNSwSypAACD+AgAYFAUAAAGRQEAYFAUAAAGRcHhtddek1NPPVXat2+vryB20EEHyQ033CA//fSTVBp1Ge//+q//kmuuuUaOO+442XvvvaVJkyb6vRkxYoT867/+q14nlx49ekhVVVXoTV0RqhzdfPPN9b72qVOn5tx2x44dcuedd8phhx0mrVq10tcMP/HEE+X555+Xcr6wVX3vV/Y2YcKEWttW8t9ZvtElNYf77rtPJk2apD/ounbtKvvtt5/89a9/lVtuuUWee+45WbBggf5HWinmzp0rw4YNM7/36tVLevbsqS/xN3v2bH178skn9XvTrFmznPvo16+f7LHHHjkfa9SovL+bqEsgHnDAATkfy9XVWX14DR8+XP+dVVdXy8EHH6y/jKgvKuo2efJkue2226TcNG/eXAYNGuR8XL0v2asfHnvssTnXqeS/s7zJoJZ33nkn06hRo0xVVVVm2rRpmd27d+v7P//880xNTY36Opw544wzKupdmz17dqZnz56Ze++9N/Pll1/Weuyxxx7LNGvWTL8vkydPrrNt9+7d9WOvvfZaptLcdNNN+rWff/75kbb7zW9+o7dT7/n//u//mvtfeOEF816/+OKLmUozY8YM/dpbtGiR+f7772s9Vsl/Z/lGUQgYPXq0/uM677zz6rxZK1eu1AVDPb5s2bJMpVD/ALdv3+58fMqUKfo92WuvvTK7du2q9Vgl/2ONUxTWr1+fadq0qd5u7ty5dR6/4YYb9GMDBgzIVJqhQ4fq1z5u3Lg6j1Xy31m+0Z6ybNq0Sf7yl7+YEaZBKgJQua4ya9YsqRRt27bV5xBcRo4cqX9u2LBBvv766xSPrPy8+OKLsn37dtl///3lhBNOqPP4xIkT9c8lS5bIqlWrpFKo8w1vvPGGXh4/fnyxD6escU7BsnTpUtm2bZvOxY866qicb9jgwYNlzpw58tZbb6X1/6jk2Sfw1En5XNQJ1bvuuku2bNkiHTt21O/juHHjpE2bNlLuli1bJuecc46sX79ev95DDz1Uzj77bH2uICj7d6Xen1y6dOlizueodXv37i2V4NFHH9Xn+Lp162a+mOVSyX9neZP3tkcD9uCDD+om6AEHHOBc5/HHH9fr7LfffqkeWym7/PLL9Xty2GGHOZv1uW777LNP5pVXXsmUe3yU66bOWV1xxRWZnTt31trmuOOO04+rSM5l2LBheh0VJVUCdV6vV69e+jVff/31Odep5L+zfCM+sqj4QwnrWZR9bOPGjfmrzA2YijGy3SpVl9Ug1Utk+vTpsmLFCt2DRr1vL730kvTv31+++eYbGTVqlN5HOVLfVK+++mr9jV7FaqpF9d5778kll1yiv/Xec889cu2119bahr/BulRstHr16tDoqJL/zvIu72WmAfvHf/xH/c1i8ODBznVeffVVvU51dXWm0qmTot26ddPvxy9+8YtI227evFmfLFXbqm++leb222/Xr71x48aZNWvWmPuz34gfeugh57bnnnuuXufCCy/MVAJ1or6+f5culf53FgcthUA/aUWd6HNR5xzCsvNK8f333+sTzJ988onU1NTIjBkzIm2v3r8pU6aYcRDfffedVJKrrrpKOnfuLDt37tTfaLP4G6xNfetX41/inmCu9L+zOCgKlnbt2tVqwueSfSy7bqX20vr5z3+uT8yrk6Uvv/yy7qEUVXYA0u7duyuqJ42iBqUdffTRennlypXmfv4Ga1O9/NTfW8uWLWXs2LGx3utK/juLg6Jg6dOnj/6pvv2qaQZyyf5RZdetNJs3b5ZTTjlF5+TqPVA9sdS0F3E0bdrULKtvzJUm+/rt1579u/roo4+c21XS32C2BTpmzJjYPYgq/e8sKoqCZcCAAfoPSEVEixYtyvmGzZ8/X/885phjpNKoE6WjR4+WefPm6blmXn31VX0yNa7ly5fX6mpZabKvX02lkjVw4ED9U01xkcvnn3+uu6Pa65Yr9TrV31rSsQmV/ncWWawzEWXstNNO8xrR/O6772YqiRrR/Hd/93f6tXft2jWzevXqxPscO3as3l/fvn0zlebf//3fTZfJRYsW1Tp536RJk3pHNPfv3z9T7rJdenv06GGmm4mjkv/O4qAoBKh/oKoPeXDuoy+++MLMfXT66adnKonqS5/9h9WxY0ddHH3ceeedmX/6p3/KfPPNN7XuV79ffPHF5kNx1qxZmXKzfPly/RqDXx7UNCBPPPFEpm3btvq1n3LKKXW2veyyy3LOfaTmO8rOffRv//ZvmXKm/t2p169eqyoOYSr576wQKAo5/OEPf9BFITtITX0ry/5jPPDAAzNff/11ppKoD7HsPyz1rW3QoEHO25IlS8x2kyZNMgO11D/wo446KnPIIYfobpjqftXqUl0zy9HSpUvNe6bmhFJ/Q0ceeWSmXbt25n7VxXLjxo05u1Eec8wxpuuzGhTYu3dvs91VV12VKXdqDqPs3059rdJK/jsrBIqCw5w5czIjR47U/6BVQejTp0/m2muvzfz444+ZSvPII484R4sGb/aEZAsXLtQzfg4cODDTuXNn/T62bNlSv5cXXXSR/uAsV+rD/pZbbtEtATX2oE2bNjoW2nffffXf1cyZM+uMZrZt27ZNf5CpDzc1K+gee+yRGTJkSObZZ5/NVNLYBPWa61PJf2eFUKX+E/1MBACgHNH7CABgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQBgUBQAAAZFAQAgWf8fT9YKrq07lfoAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax1 = plt.subplots(ncols=1, figsize=(8, 8))\n", "ax1.imshow(data.data[50], cmap=\"Greys\")\n", "\n", "im1 = ax1.scatter(\n", " features.where(mask)[\"hdim_2\"], features.where(mask)[\"hdim_1\"], color=\"red\"\n", ")" ] }, { "cell_type": "markdown", "id": "163a7045", "metadata": {}, "source": [ "If the features have the same threshold, tobac keeps the feature with the larger area. Otherwise the feature with the higher treshold is kept." ] } ], "metadata": { "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.13.11" } }, "nbformat": 4, "nbformat_minor": 5 }