Spineweb_threshold.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nibabel as nib\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "from torch.utils.data import Dataset\n",
    "from torch.utils.data import DataLoader\n",
    "from torchvision.transforms import ToTensor\n",
    "import torchmetrics.functional as evafunc\n",
    "import torchvision.transforms as transforms\n",
    "import numpy as np\n",
    "import torchvision.models as models\n",
    "from segmentation_models_pytorch.losses import DiceLoss, FocalLoss\n",
    "import segmentation_models_pytorch as smp\n",
    "import os\n",
    "import imp\n",
    "import CTimageLIB\n",
    "imp.reload(CTimageLIB)\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 自定义数据集类\n",
    "class SpineWeb15(Dataset):\n",
    "    def __init__(self, data_path, transform=None):\n",
    "        self.data_path = data_path\n",
    "        self.transform = transform\n",
    "        \n",
    "    def __len__(self):\n",
    "        file_list = os.listdir(self.data_path)  # 列出文件夹中的所有文件和文件夹\n",
    "        file_count = len(file_list)  # 获取文件数量\n",
    "        return file_count\n",
    "    \n",
    "    def __getitem__(self, index):\n",
    "        metadata = np.load(self.data_path+str(index)+\".npy\")\n",
    "        image = metadata[0]\n",
    "        label = metadata[1]\n",
    "\n",
    "        image_data = torch.from_numpy(image).float()\n",
    "        label_data = torch.from_numpy(label).float()\n",
    "        \n",
    "        if self.transform:\n",
    "            image_data = self.transform(image_data)\n",
    "            label_data = self.transform(label_data)\n",
    "        return image_data, label_data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_pixel_accuracy(prediction, target):\n",
    "    # 将预测结果和目标标签进行逐像素比较，并计算正确预测像素的比例\n",
    "    correct_pixels = (prediction == target).sum()\n",
    "    total_pixels = target.size\n",
    "    return correct_pixels, total_pixels\n",
    "\n",
    "def calculate_iou(pred, target):\n",
    "    intersection = np.logical_and(pred, target)\n",
    "    union = np.logical_or(pred, target)\n",
    "    iou_score = np.sum(intersection) / np.sum(union)\n",
    "    return iou_score\n",
    "\n",
    "def calculate_dice(pred, target):\n",
    "    intersection = np.logical_and(pred, target)\n",
    "    dice_score = 2 * np.sum(intersection) / (np.sum(pred) + np.sum(target))\n",
    "    return dice_score\n",
    "\n",
    "def visualize_results(image,label,output):\n",
    "    fig, axes = plt.subplots(3, 1)\n",
    "    axes[0].imshow(image, cmap='gray')\n",
    "    axes[0].set_title('image')\n",
    "    axes[1].imshow(label, cmap='gray')\n",
    "    axes[1].set_title('label')\n",
    "    axes[2].imshow(output, cmap='gray')\n",
    "    axes[2].set_title('output')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of the dataset: 365\n",
      "----------------  Loading the dataset and evaluate  ----------------\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJYAAAHWCAYAAACVLnA6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABTaUlEQVR4nO19eXhU1d3/Z5bMlkkmCZANEvZNtNA3QKTFBU1F27pUH1Hq2wfRSrVgrUitvK8S8KXFuqBVEVFbqFpxeR+Xvq1LLYuiAiqgVBAEDCECIWSbTJaZycyc3x/5fU++984kTCA3MyPn8zz3SebOnXPPvfdzv9v5nu8xCSEEFBR6GeZEd0Dh2wlFLAVDoIilYAgUsRQMgSKWgiFQxFIwBIpYCoZAEUvBEChiKRgCRawYWLNmDUwmEw4ePJjorqQsFLEUDIFJjRVGIxwOo729HXa7HSaTKdHdSUkoYikYAqUKY0BvYw0ZMgQ//vGPsXHjRkycOBFOpxNnnXUWNm7cCAB49dVXcdZZZ8HhcKCkpAQ7duzQtLdz505cf/31GDZsGBwOB/Lz83HDDTegrq4u6tx0DofDgeHDh2PVqlVYvHhxTMn5/PPPo6SkBE6nEzk5Obj22mtRVVXV6/fjpCAUorB69WoBQFRUVAghhBg8eLAYPXq0KCgoEIsXLxYPP/ywGDhwoHC73eL5558XxcXF4r777hP33Xef8Hg8YsSIESIcDsv2HnzwQXHOOeeIe++9Vzz11FPitttuE06nU0yePFlEIhF53Pbt24XdbhdDhgwR9913n/jd734nCgsLxfjx44X+US1dulSYTCZxzTXXiCeeeEIsWbJE9O/fXwwZMkQ0NDT0xW3qFopYMRCLWADERx99JI955513BADhdDpFZWWl3L9q1SoBQGzYsEHua21tjTrH2rVrBQDx/vvvy32XXnqpcLlc4vDhw3Lfvn37hNVq1RDr4MGDwmKxiN/97neaNv/9738Lq9UatT8RUKowTpxxxhmYMmWK/FxaWgoAuOCCC1BcXBy1/+uvv5b7nE6n/N/v96O2thZnn302AGD79u0AOhyGf/3rX7jiiitQWFgojx8xYgQuueQSTV9effVVRCIRzJgxA7W1tXLLz8/HyJEjsWHDht667JOGNdEdSBVw8gCAx+MBABQVFcXc39DQIPfV19djyZIlePHFF1FTU6M53uv1AgBqamrQ1taGESNGRJ1bv2/fvn0QQmDkyJEx+5qWlhbPJRkKRaw4YbFYerRfMGd7xowZ+Oijj/Cb3/wGEyZMgNvtRiQSwcUXX4xIJNLjvkQiEZhMJrz11lsxz+92u3vcZm9DEctgNDQ0YN26dViyZAkWLVok9+/bt09zXG5uLhwOB/bv3x/Vhn7f8OHDIYTA0KFDMWrUKGM6fopQNpbBIIkidOHCRx55JOq4srIyvP766zhy5Ijcv3//frz11luaY6+88kpYLBYsWbIkql0hRMwwRl9DSSyDkZmZiXPPPRf3338/2tvbMXDgQPzzn/9ERUVF1LGLFy/GP//5T3z/+9/HLbfcgnA4jMcffxxnnnkmPvvsM3nc8OHDsXTpUixcuBAHDx7EFVdcgYyMDFRUVOC1117DnDlzsGDBgj68ymgoYvUBXnjhBdx6661YsWIFhBC46KKL8NZbb2m8PwAoKSnBW2+9hQULFuCee+5BUVER7r33Xnz55ZfYs2eP5ti77roLo0aNwsMPP4wlS5YA6HAkLrroIlx22WV9dm1dQQ3ppACuuOIK7Nq1K8ouS2YoGyvJ0NbWpvm8b98+vPnmmzj//PMT06GThJJYSYaCggI5rlhZWYmVK1ciEAhgx44dXcatkhHKxkoyXHzxxVi7di2qq6tht9sxZcoU/P73v08pUgFI7CD0448/LgYPHizsdruYPHmy2Lp1ayK7o9CLSJiN9dJLL2H+/PkoLy/H9u3bMX78eEyfPj1qyEMhNZEwG6u0tBSTJk3C448/DqBjmKKoqAi33nor7rrrrkR0SaEXkRAbKxgMYtu2bVi4cKHcZzabUVZWhs2bN5/w95FIBEeOHEFGRoZKHe5DCCHg8/lQWFgIs7l7ZZcQYtXW1iIcDiMvL0+zPy8vLyoQCACBQACBQEB+Pnz4MM444wzD+6kQG1VVVRg0aFC3x6SEV7hs2TIZXeaYMWMGjhw5gi+++AJCCFgsFoTDYdhsNrjdbowdOxajRo2C3+9He3u7zAowmUwwm80ys8BqtSItLQ0Wi0Ue197eDrPZLNskyWgymeBwOOB2u2Gz2WA2mxEOh+WYnclkQjgclr8Jh8OwWCxIS0tDOByG2WyGEAKhUAhtbW0QQsBqtcoxRToP9ZXaMZlMCAQCCIfDiEQi8jd0fDgcht/vh8lkQlpaGhwOh+xXWlqaPE4IIdumz62trQCAUCgk763JZJL9tlgs8p6sWLECGRkZJ3xmCSFW//79YbFYcOzYMc3+Y8eOIT8/P+r4hQsXYv78+fJzU1MTioqKkJaWhmAwKG+UEALhcBhWqxV2ux0DBgyAzWaTN1gIIckCdDw8s9ksyUIPMBKJyIdmsVgQCoU0hMzIyEB6ejqAjtwnIQTa29sBQP6eCBWJRBCJRJCWlgaz2awhtMPhQEtLCyKRCGw2m4a8PJ2G9z0cDiMUCsmHbrVaEYlEEAgEotq32+0Ih8NwOByS+NQW3QMibiQSkfeS7gltNpsNbW1taGxslP07ERLiFdpsNpSUlGDdunVyXyQSwbp16zRZmgS73Y7MzEzNRr/x+/3yhtMbbbVa4XA4ZNIdEY9uGIGTh95kLsX0uU4kKex2O6zWjneSHhQ98EAggFAoJMlGU8mIeHQ8PTC32w2r1aqRFgQ6jtriv6P91G/ql81mg91uRygUkveGrslisUjymc1mWK1WSR76nsDbJ9Jyc+RESJgqnD9/PmbNmoWJEydi8uTJeOSRR9DS0oLZs2fH3UY4HEZra6uUEKQGTSYT0tPTYbPZ5I3n4DeNbjapQUB7g7n6FELA4XDIB0KEpgdDhKLz8d+1t7fDarXKB0376SXw+/0aKacH9ZdLLyI0nYu/IGazWRKdq2J+7kgkoukvSTWr1SrvJ90Ls9kMu90e97NJGLGuueYaHD9+HIsWLUJ1dTUmTJiAt99+O8qg7w4kHYDOt9tut8NsNiMzM1MjCUjcc4lFKorUJQD5JvMHSA+AjqW2uPQhYvIHRvYZkcBkMiEYDCIUCkn7jOw7UkUkVYkcdH6+n1Q1EZUIQ+cDIFUktc1tNbpf3C6kvySJ+f90XE8yUxNqvM+bNw/z5s076d9ztUMPwWazAQBycnKkwUngNxWARjWS0U4PXC8NaKOHyaUHEYlccPod2Wb0UIPBoFSPRDySfiRtiHQEThTqCzkWJFWATrVOkpT6k5aWFmWU01/+UlAfuXFP96m9vR0mkwkulyvuZ5PS2Q3t7e3yptHbTwRzu91RRiZ/+ECnyqM26EHzY8k4pwfA2yJwG476wo+hh082oMlk0qhNIiw3orlnCXQSglQmnzBBLxARlnuWtJ/a4oY5tUu/oWvn56P/6eWJFykRbugK5K4TyDZwOBxwOByIRCLyAXKpwyUWqQyr1aqRRtQWGfH0G25zkYqh9unmc8kBdO9FcQ/UarXC7/cjFArBbDbLayBVV1NTg+rqaqSlpWHIkCHSM6XzkpFNpON9p2vhUpuulT5z04F7znSvehKMTmliBQIBjSdHDzUzM1PaFlyacKOYJBDt5zed/vLQBD0YOg+XRpxc9Bsimt/vlwY5fzhcqnDvjB4u3xcOh9HS0oLPP/8chw4dQr9+/RAIBDBs2DBkZGRIm5BeBlLBRKy0tDQp3bmko7atVqu0D0ktkrqldjkJ40FKE6u1tVUT8yHjNzs7G1arNco95jeHk41UYCzvMRgMauyw9vZ26W0CnQFH7tKTQc5DDiS9KC5G5OMOAIAoFcfVWCQSgc/nQ1tbm1SrI0aMwIABAzTB3EAgoFFvtJ/Oy++F3ingtiX1l/rCVeWJkNI2lt/vB6CNdrtcLuTl5WnUFIGHAehBkyShv4DWuCdwG4fICCDqLz0YPWEIFBLhxjg/BzkgvH90vv79+yM9PR1msxlerxd79+7FN998I6WyxWLRvCREJn4tPJCrJzSFKEhykbSi77gjdCKktMTy+/2w2WzS8wmHw8jIyJABVO6GE3nozeOGtl5V6tUVkYCMYVITJBGoXToft+14qIH/paESLvm4F8m/ozYGDBiAoUOHorGxUcaUKFrOpSOdm9tLFCYhD4+ulTsI3Oak6+ehDy7JToSUJ5bdbpfESUtLk8M4dMP5GBzZXVwFkBfJbzCBVJje6G1tbZXkpHaJFHQuiiFZLBYEg0G5H4BmSIYeIhFSLzV5fzMyMtC/f384HA7k5ORIKUU2GUkcTgSSOgCkrRcMBjXjhzzGR9dCKp6/nD2ZtZ3SxKLYEz18l8uF3NxczU3lRrLeKwSix+U46OZyiUcIBAKS1Hq7SG+oOxwOhEKhKDXDJSW3f6gdIgrZT2lpafB4PHC5XPB4PFIqhUIhSRZ+TXSt1AaBJBl/OUiFc4nGg8Y98QiBFCcWvU00gOx0OpGVlQVAO1LPpROXBPQmc9Jw26arm8nHJqmSDEkf8ur4b3nglnuO9DkcDsvYFQAN6WggmdtaNpsNNptNI5koNsYNbJJmXLVx75agv1Y+LETk5oPk8eBbQSxSOW63Oyo6zKUO3wdojXMAUapIT0a9MU4DvURsHkjVqxce1eeRdiEEgsGgHNPTB2JtNpu06xwOB7KysqSNR+qfpAv1lSQsXRP1hcISPN1GL631DgsPjOrvV3dIaWJx45IyBShmw4/hopy74dygp2P1MTF97AvQepLk9vO0G6BTCvLALI83EeHb2trQ3t4uPToiFhnh9NLoh2/IeOdEpHPa7XaN1KG2eaSfR+CJOCThiLSkHomop43E0t+kEw3j8Icey7aim9xVaSL9uan9UCiE5uZmpKWlyU2v1jhhgA5vrrW1VaO+OPkBaGwxIonD4ZCSivpL3ip5m7E8XyKYfiRCH+wlE4LH5boqbNIdUppYXPKYzWa4XC7NAHEsu4qI1hV59O61nhB6dcg/+/1+mcVJ+VF6SUlSimJwFHeKNcTCpYndbpeEIcnGpRiPsnNVzPcD0IRm+EtH++jc3CnqSWCUkNLEogsnklBEnD9s/nB4wI+3AWglEH8z9f/TA41FLtooR5+Cqfw4HkUno56fU99fHkvTh0+4iuWRflKtJIHS09M19hWpf5L45F2Toc81AYGr+bieTdxHJiH420RvJx+aADoDgPqhHA66YbFEPhGJQG3rh0PoWJ4KQ94eDSzz4CUdrycbtc2TCXmfeCxML63C4TDa2to02RH8Wrh3TH2J1S+6F3TungRG5bPp8S+SCFyM0zAJf5MJXHpxj0jfFk8t4YlxesnB40J6D4+PPdLAL1fPQGdGg95m4UQi6aHvD8WagM4ALj14rrqBDglO6clcJZMdRcTiMS+6RvJg+XDXaeMVcs8lVkiBbjK/OaQGuUcYS4LRA6UxMjoXj/lwT4uThtrj42z0nT78wUmq97x4Sg9lqJKXxtOMSTLxNB8e4KRr0qvlWJ9JPZNq7Um0nSOlicW9Ir2qIMRSgfoxOiD2YLS+Pb13CUAj2ag9bnOR9xYr8Mr7zgO0emlI10kvRawIPffiqJ1gMKh5Gahv+n10Lh4w5iYA91LjRUoTC+gcQtG76vQ/9wS5vaKPXelFPZcGfCCZjtVLSsqFoja5V6pP/dWra1JnfAID9xC5iiYjm3+vl2wAZMoQdzY4YTix+HVzSaonUk9srZQ23knc6wkUy3bh5OpKvOu9pVhRaU4ObmzznCt+vB7Uv67UL9lPPLjK+8ePBaCxK/mwC49jkR2l7yN3PnhIQk8gLtXjRUpLLLpZ3BOM5bYDWjUWayhHr/K4KujOGSBpQ1KDqxzeLrXHE+dijVFyR4DCBvpYnD7wSUSN9ZLRvaFB6ljjgtwW5W1TW3zYKF6kPLGATtXGp3rRX/0b2tXN4SSMlUJD0HuHAKSq4+EKTjBOGt4Gfxm4Q0ARcn2Snv569H2mNqk/XILrbTlSdaQ+9UNc+mEi/flOhJQmFr9h/Kb15AbwG8dVTazzEHnowZPaIgNdPzZI6I6odD5uk1FfSGrpwwXUbz159Y4Hlz704nGHgNtl3M6kc3GyxlLf3d7XuI9MQtDbRf9TQp3eSI9FNP0sFL1K4vt4MJRyn/x+v4xk6z1M+p/vj3UMP5Yb+iQxSH1yKcvtRB5n48NU9D+3E/1+v5wsS8fw6+ESjv+O9ncV/+sKKU8soNN+oaopQGwPRi9FCLHIpP8df/D0GzJ29XYLDzrSuXiuOYc+hYXbdXxyBEXv9c5FrPAFPy+ROhgMSgJRrIqum4aB+AvAr727l6IrpDSx+NsViURkKWu9SNffkK7UBf/LRT8PC5B04UNJfJyNG8P0W95Hgj42xF17TmBSuyRtiKScgLw9ciK4sa/36ng/eIiiK6+P+h+rcElX6HViLf7/y8zybcyYMfJ7v9+PuXPnol+/fnC73bjqqquiyhnFC/7GCiHQ2NiomVrFwSPLnCixHjagnXfI3XYetqDhEfpLUoFnD1AfMzIyNEU19FKAvFs6N1dBFP3ndhCXRvq4HHdeKBQCdEpHSu3xeDzIyclBTk6OzHrgVXG4SuZkjgeGGO/jxo3Dv/71r86TsLSL22+/Hf/4xz/wyiuvwOPxYN68ebjyyivx4Ycf9vg83MYAgJaWFjnvL9ax9GbGijnFioDznG9SU5SlylUvJzMnOg+HUBu0n9x/fXoN0GlzEUHN5o6UIJOpI3NCH2MjkuuNeBoGam9vR1pamiyYQi8XEYwcEp7KQ/3mA9UJj2NZrdaYBdS8Xi/+9Kc/4YUXXsAFF1wAAFi9ejXGjh2LLVu2yFVHTwZCCLS0tKC1tVXWr6IINdlHnDixgp9dgQhCG2WN0iAzxbFsNpv0vPSShFcKjEQi8Hg8KCwshMvlkiktlOoMQCYPtra2amo+6J0Gfj28zgRXjZT2zIdsOHlIyvL4HZdSPBQSLwwh1r59+1BYWAiHw4EpU6Zg2bJlKC4uxrZt29De3o6ysjJ57JgxY1BcXIzNmzd3SSx9DdKmpib5P7+JgUAAzc3NyMnJiRqb03t9emLpVSKRhXua9JkMaZolTeEGm80Gl8slE/4I3J232WwYNmwYBg4cCL/fj5aWFkkcnt1ps9lQUFAAADh+/Djq6+s1g+H0l/rHp3/RNZAUJWkthNAMPYXDYTQ1NWmmiQkhZB4/T3Ak6Rovet3GKi0txZo1a/D2229j5cqVqKiowDnnnAOfz4fq6mrYbDY5k4aQl5eH6urqLttctmwZPB6P3Phyudx4DQQCaGxsjHpre+LNUJv05pLq4HnjXJ0QUa1WKwYOHIihQ4fC5XLJWlo2m03aSRkZGRg/fjyKi4vR0NAAr9crSUgEIYkWiURkbYqCggIMGTIETqczao4jAElySnXms3YASKkIdHrDPChLxURoRhCBVCVlr/YEvS6x+MLY3/nOd1BaWorBgwfj5Zdf1iy63RN0VYMU0NZgCIfDqKurk3Pj9MM6fLiHg2wfPljMVSN/w4lsFosFTqcTkUgEWVlZGDFiBM4880wEAgHU1dVJ+4vbO3l5eRg4cCAOHTqE1tZW+P1+OfmDG+r0O26f5eTkAOioWEw5VlxqkqThbQCdXiIfOOcOAL9/JBGBzqn4fIA74cY7R1ZWFkaNGoX9+/fjBz/4AYLBIBobGzVSq6uitgS73R6zTCG3ZYCOB19fXw+/3y9rHPTkTeMxIP3veJiBpEt6ejomTJiA0aNHIz09HdnZ2WhubpbFZHngMhLpmD3T2NgoJRHPeCUJQ33m9gwROicnB62traivr5cEodQYmrBLpCGjndtHXJKT+uPzDun31CZJv5NZvNzwOFZzczMOHDiAgoIClJSUIC0tTVPUdu/evTh06FDMorYngj5eRedrbm7WxIKA6Og3D3QCsWtY8YkL9GZz7zItLQ0jR47E6NGj4XK54PV6kZmZif79+8NkMmnCD2azGT6fT9pePF4WS7qSjccHup1OJwYNGgSHwyFLbrtcLjgcDqSnpyMjIwNutxvp6enweDzIyMjQeMi8/UAgIPtGpKeQA0lwUuX6Iaq4nk2Pjo4DCxYswHvvvYeDBw/io48+wk9+8hNYLBbMnDkTHo8HN954I+bPn48NGzZg27ZtmD17NqZMmXLSHiFXIeQZ6stG6x8aj3DHyhagduk4/fgdSaKWlhZUVFSgra0NXq8XDQ0NCAaDyM/Pj4rMB4NB1NTUyIp8XFrwuBodz/sPdHpkWVlZKCoqkjaT0+lERkYGMjIy5EtAjobT6YTT6dQY/Bzk9ZEdpbfPuNeYcBvrm2++wcyZM1FXV4cBAwZg6tSp2LJlCwYMGAAAePjhh2E2m3HVVVchEAhg+vTpeOKJJ07qXDxASMQJBoOor6+X7jefus4nhuonGhB49Lyrm0mqKhAIoLq6Gl6vFzU1NWhqakJ2djbq6uqiYldCdBQTqa6uRnZ2tkbtkUenLx6it/vItsrKykJdXR0aGxulB8kzK/h18aJpXLrHmmNInjfFv0hN81BJvOh1Yr344ovdfu9wOLBixQqsWLHilM/FHxrQGSaora2Vs4j1DxfoNMbJ9gFiR+H19pY+w5MeqslkQnNzMwDg6NGj2L17d9SoABGppaUF2dnZUTEoPrhM/aI+k4cKQF5Xbm4uvF6vJpzCJS4fQOdGOycy2Xhkv1JA1GQySU+QpzefNmOFsQZNI5EIamtr0dzcrAkK0nfcs+FE4v9zlQJ0DreQBCEVl56eLss2BgIBOBwO7Ny5U4Y8qI9Ap/caCATg8/mkdKKHRgY9JweXXDTc0t7ejkAgIG0q6h8/H5GND21xqUUko5AIdyICgYD0UqkP1LfTJrtBPwxC8Pv9qKur0xzLxTpXQxxcalD71C4FN8m7slgsOPPMMzFhwgTU1dXB5XKhvb1drrdIkkafwkJhAb3Nx//ng8Z8kJt7agCkpOHqjPoKQJKRf88JR//TC8dfJjqn3gaNFymf6MffLPq/vb0dDQ0N0g3nLjUfatF7ejyoyNslyUdtWCwWjBw5EpMmTYLX60V9fT1cLhcqKirkFHtCrMkYNLaoj6Tza+C2FbXDnRGgczyTfuN0OjU2Z3t7u8z44OoS0DoI3J4iCcbrbjkcDgAdY7Hx4ltDLL4vEomgubk5KiBK3/FiZrEGnek3FLeicUC73Q6Px4OCggI4nU589dVXaG1thdncUVfe5/NppAwnLre36Lyx3HgeeuDf8cFl7q0R6UOhEOrr6zVxLCr263Q6o8Yv9RNmaeyTv1x0j/TB23iQ8sTi4A+EVoGgspH8N7EeGFcFaWlpcLlccLvdyM7OljEhIQS8Xi/q6uoQDAZlDIkMeqrNxb0vbpsIIWQlGq5igejJrxTHIonW3t4uKwhSX+lYnrJDv+V90g9xUb/onnFVxx0VatPtdkcFbU+ElCYWB5dO9DbrQwrcI9KTkmYtp6WlydiQ2+1GZmam/AsAHo8HxcXFaG9vR3NzM5qamhAKhdDU1IR+/fqhoqJCk+FA4EQj0vDvSRIRcfQqUm8889KTNICtr2LIqzNzFav3SDlI4tELlJWVFTU2GQ9SmlixLlj/MPUPGNCWieS2F3ldfr8f9fX18ns+A5mi2gUFBSgqKkJhYSEOHjyIY8eOobi4GIMGDUJFRYVUKVwF8Yg9DQ3RGBwfm+NhAfqdvkoOZVhQ30lNcvLx6f38ZeKD7GRfcTPAYumojtjY2KhxaHriFaY0sbgRzg1bso34wCugrbaitxf4m6svpc2/a21txfHjx/Hll18iNzcXkyZNQl5eHo4fP46mpiaMHTsWDQ0NMngZS2Lqh0moXzztmM+e4TYRbUQsUnN8EgZ5g9QGt8+AzngceZg+nw9CCBmlD4c7atEPHjxYqm+yTeNFShOLgw+1mEwmTQIdf+PiFef6WBIPMNJNPnLkCDZt2oTp06ejuLgYFRUVsFgsmDhxIrZs2QKfz6cZ9KXxOKrKR1KNp7RQ4JKXy7ZYLDLUQXEsXhqJe6uxVnLlkl2IjmGvyspKucoF0LFknMfjQX5+viS53rk4rVShPtZE/3s8HulK6411fqO51DjREA4/jqLTLS0taGhowNixY2V8avjw4fjmm2+wZ8+eqGClxWKRqcik5sh7o/NwO4i8UQp3UN4VLalHfaa6oxRi0Re9pev1er3Ytm0bdu3aJRe56tevHwoLC9HS0gKr1Yr+/fvLa+fB4tNGFQLRbyPQYafwtBw+o4UnytEN44TSx3loH1c5/NwWi0UW9B89ejT8fr9MSCRy8GETyhggScWDl7HUM1ffpMIol4uug48u8BeFIv38Hh0/fhz79u1Dc3OzJF5DQwPM5o68ervdjoyMDBlqICeA0p7jRcoTC9DOeCE1SKuBcrtAP1QDRE+i4IPQHPpxQ7PZLMswUgJjv379EA6H4fV60dTUpLGPyEg3mzsK1Pp8PtkukUAfhiCVxAulUeYp9YPSYniaDY+j6WtK0F86BugIKFO6N5kOJDkpRQfQEvhESGli0cPlN81k6lgPmmJK+gfFA6P66DaBJFGsNF46J9DxoAOBgFQrNpsN6enpCIVCclkUQFvog9RNenq6NPABrbolo5xXjwmFQqirq5PhDT5/kJOF95dIwgOeLpcLLpcLdXV18n5QunhxcTGKi4ulGuVpySQt40VKEyuWQW4ydaxXSIFLMm75BFKg6wIh9FC4J8nnCOqPbW9vx549e7Bnzx5YLBYZUCXDnaQjpRC3t7ejsrISI0aMkAsHcEeBYllkiFO/6+vr0djYKPO5SMqdaEhIT9hgMCilHKnovLw8jBkzBkOGDIHb7ZaeIBGMBqZTYhX73gS3Ici+4uQgScbXmtEThw/vUFs8us2HN4ho+vThSKRjFdSampooz4yrX5o4MnDgQKSnp6OtrU2eh8eaKJ+ssbERDQ0NUrryibQkTfQzr3m+PEk4ygej3HibzYaBAweipKQERUVFcrIGoF21lWJdp03knT94+myxWODxeDQqQO8V6mNL1A7fpyeG3sjmA71cbXGpwe03evgUYqiqqkJ9fT3y8/PhdDqlV8fXcCZbjUIC3AbjXiXvG3+Z9M6A1+vFgQMH4Pf7MXDgQBQXF2PUqFFy4VAeP+NDQ4TTZhCawO0Jp9OpWa8QiF19BdDaItxT5ODEpXPQfv3D1A88E4igfJxQCAGfz4fW1lY4HA5NvXcKcpLK47ORuSerj6pzNc/7TvYYpVBnZGTgP/7jPzBkyBA58YNCF9zYb29vh8/nk4Q7bRbCjCUhnE6nvFmAtpRRLCOdPyR9QJH+59KJjF3aT+Nx3BEge4r6Q32jMAMPaNJ3FMsiFUYZGEQqMuhjlQmIJZH5GCMASZL8/HyMGDECw4YN08wjpKlofFo+j963tLRIyRkPvhXE4oRwOp3y7ecShG643paKBR5hJ/XK33xORE5anqnAPVC91NGvnQNoh6e4U0JqlMjNI/+8tHis8VCuKoUQcLlcKC4uRmFhoeY+8PRjPgBObZL3e1pJLK4GKIbFiaV3w2MRAdDGp7hE4SEHPvmTwLMxY9lg1HasgdxYUpR/Fw6HNV4j/w23tbhU5dfDJ2lYLBaMGTNGhkF4lkWsWBu1QY5FrAWwukNKE4vbOXQjKWeJSzH9G83tqe5iM7F+o4/Mc4Lq7S9Oaj6Lhs7Lh1r4NZFKpGlXPAZHJOBjjkBn8TRO3vb2drS2tqK2tlauPpuVlaWp4kfE5eW79USmEYKeTFxNaWIBnQY1PRw+Dqc/Tu8yc+OeR965mqPvueSi4zhBueHLVSIdqz8PD9RSG9yW49X3uNNAcwWpaAd3XEgyAR020/79+/HVV18hEolg/PjxckSCpBWpOLpnlBhJfSKvUB9ojQcpTSz+kMgWogFbPtNELxl4Mp1+TJDbTdw+imX4c/efgpqx2uPnMpvNcn4jJy+FIfTlIPl4oMvlgtPplP0islO2LKmzYDCIo0ePYuvWraivr8e4ceOQk5Mj16amzAgiCqlEmpNJqcyUPgN0DnLHi5Qmlt6To7ExTiQuzbitxT/HGifkTgG3Z3iaCldRpC64g8BDFVwScknGY0XcSOekAiDToPl5KIBJ0/bJGzabzaipqYHP55OD4xTbq6+vR01NjWyPT1wlolGhNlKR/Jh48a0hFqAljN4wpjgS/a//LYHbMnwff+CcUERAUkOxbj5PtiPCkJ3D7TEiJH8ZKAJPkzp4nIuyXUmyUP94puuZZ56JwYMHSyne1NQkB8ipwC2di87Bs1OJWHSN8SKlicVXWuC2DfeU9J5aLPuJpAgfE9QP6XBbhgjLwxlEPL3hy/tF077IaOagNrjtRueipD/qD6XOUFyNshqon1arVdZ9zc3Nhd1ulyEFyumnNawBaIKffBySxhVpGhllZMSDHk9Yff/993HppZeisLAQJpMJr7/+uuZ7IQQWLVokp0iVlZVh3759mmPq6+tx3XXXITMzE1lZWbjxxhvlFPWegBvi/MHoHxCPAelTU2KpO04evX0Wywuk/XQcl3Y845MePqk/fbiD9vG2SAKRGqZJHEQqPtuZS0u73S6Ha3huO9AphbjnyQOv1EduZ3EPNB70mFgtLS0YP358l7UX7r//fjz66KN48sknsXXrVqSnp2P69Oma0onXXXcddu3ahXfffRd///vf8f7772POnDk97YrGVqGbwfPV9cM0nBxdQX/z9DYcP58e+nbpgfIJC7x/QGepJB5GINXHidje3o6WlhaZTsyDr3RtkUhEepSk6oDOor8k/QKBgOwDHU99oZx8KpGUlZUFu90uryNe9FgVXnLJJZqqfRxCCDzyyCO4++67cfnllwMAnn32WeTl5eH111/Htddeiy+//BJvv/02PvnkE0ycOBEA8Nhjj+GHP/whHnzwQRQWFsbdFxpO4fGbWC4xV2V64xmIVoexoB/6iUVO/XgdryxI0XZ6OGTf6CP2egOf2m1tbZUE41PrucTlM5rpeqm+KZkNdrsdNpsNzc3NsgI0xbKoH3QPiJjNzc1obGzskVbp1doNFRUVqK6u1hSv9Xg8KC0txebNmwEAmzdvRlZWliQVAJSVlcFsNmPr1q0x2w0EAtLopA3QTh6gB0hvo368T6/S9JFq/Xd8n94JALS1tPh56BjuPer7Qv0myUAJdfq0Z6AzLkfGPndQSHLyfWTUU64XgafcOBwOtLa2yrE/m80m1SJ5h5RlIUTHjCCXy6XJhT8RepVYlGeUl5en2c+L11ZXVyM3N1fzvdVqRU5OTpcFbrsqbktvOP2lIJ/e24vlOXKjXh974jYYJ4u+TSKQvm2gczFwMq4pDkVSgiQHhUd4KSFOKopb0fXyVcd4IJgPXvOgpt5TNps7UqMjkYickEGBUj4ITsY7VbZxu909qiGbEtVmFi5cCK/XK7eqqioAnTeWapSS+qHEtlg3lq+8oCeVXsJwR0Dv7QkhZGASgIbg9IC5reJ2u+X/6enpcDgcMo7EpSkfZqL/w+GOMpE5OTnIzMyU0XEePecShvrA+07X3tjYCJ/PJ5MHyfalNimEwSfI2u32qJLoJ0KvhhuoQO2xY8dkjXL6PGHCBHkMlfohUEGLrgrcdlXclmIufBaJ3++XdggfGNaP3/EAKBAdbOXQe4882k62CwfNcKaMTDLGYxGU+k2pNmSPUU0qqr/Ao+0Oh0POsiESU7/1A+c8EZBeglAohJaWFo207devHwYMGCCNd273UcZpenp6zOcTC70qsYYOHYr8/HxN8dqmpiZs3bpVFq+dMmUKGhsbsW3bNnnM+vXrEYlEUFpa2qPzUcGLtrY2tLa2ataDiZURQOAJcwS9lOIBUUBrU+kDpPq2zebOqVTcdtJ7qZzoRDIatnE4HLK2KE1WpXNS+3z8LlYogDxKmuRhs9mQmZmJwsJC9OvXD+3t7fj666+xZ88e7N27Fz6fT9N3UuF2u132K170WGI1Nzdj//798nNFRQU+++wz5OTkoLi4GL/+9a+xdOlSjBw5EkOHDsU999yDwsJCXHHFFQCAsWPH4uKLL8ZNN92EJ598Eu3t7Zg3bx6uvfbaHnmEAKT9QWqA3HaKGVEaCA9kAideSDuWxNLPDNYb5ECnl0YGOXlVPIDLCUvHU//4AgH6VGe9NKWwBKlcmhlEZAI6xy958NZs7qgZP2TIEBw5cgT19fUQQqCpqQlHjhxBenq6xnblHqqhg9Cffvoppk2bJj9TYf9Zs2ZhzZo1uPPOO9HS0oI5c+agsbERU6dOxdtvvy2LdwHAX//6V8ybNw8XXnghzOaOQrePPvpoT7siaw0QSC1QzSo9QfgD1QcoY2U48Hb1IQD9b7iKpTE92h+rD0An4SKRiJQOPKxBx+tzz0l60WQLMur5DB4iHrVBtbvIbKBBaSojTl4n1cjns4UIhhLr/PPP73Yw0mQy4d5778W9997b5TE5OTl44YUXenrqKDidTk2NA7o5LS0tUdF1PhZH4GN3entLj1j7uEolw5mcCDoXT57T23x0Tp7qzLMpKMxARjWPsgOdpZfItuLXSlKZ+kezdGh/WlqaVLO8pCSRlUIMJOW4IxQPUnqskE9C4DZSU1OTNFS5d0RE4C465XnrA6P6KD1Xpdxzo88k0Ugi0H4+nETtckOaBpQJPMWG1DipepJkRCIiInmHvH4DOTJ8vWohhCz6m5mZKSUdTaAVQkipSfMyycOm+xwvUppYfAyLx218Ph/a29vhcDg0dd4JpF6IAN2RittXeu+PSMffdm6o6ye68qAp90ypL5xQ9HuSbHqpy6PkNCSjDzfozQGTySTnXBJh0tLS0NbWJivh0Dn5TCE6Px+WOxFSmlhpaWlyzRwhhBTrDQ0N8Pl80hjmcSH6y4nFHxRXJ10NuvJEQbJ/KC6lzzQgovCNzqMPvnI1yUsbcYNfTxR94TQ6Ny9+IoSQ/aPPQMeIBq2okZeXh/T0dLkMij7bIRQKwev1xv1sUiJA2hWsVqucEk4qIhwOo6WlRQ778IehT1bjD5lD/xC5+tQHXXnBfT4GyaUb2Ul6h4GCm6QKOUH4DGQu3Xj4g6QJqSySdPyFAbRZIFyykg2WkZGhiVFxx4H397SZYm8ymeB2uzUJdqTmjh49imHDhslgIyGWq08qRj/ArB+u0ROHxtYyMzOjFirgao2rQ/7XYrFIr47+5xKOtwVA2jpk7/B0GH2SIX8ZKHxB/5Oai0QimpJFRFCa2Go2m+HxeKTtRsvWxIOUJxZ5Nn6/X2Mb1NTUoK2tDenp6ZrxNT6IzD0kQOvl6Y1tbgCbTJ1VmUlS0cPXE5OkJCcU2TIUjefRd+oTJ5Re8lB7NODOp/nTd3TNlKxHcT0CjQaQ1PT7/QgEAjh8+DC+/vpruQB8bm4uxo8fj379+p0+xjsAWfeA5r9xA97r9SI9PT0qnsUfPhmxXPXpSUUgKcHVEnlW+tAFt7NigdQMEZMeMJc8+hAInY9+T0Y7D0PQy0VSiOYD6kcahOgoqVRbW4vGxka0tbWhpaUFR44cQW1trVwQiqoHTpw4UZYuiAcpTSx6S6nCL48XBQIBHD9+HPn5+bLGE0ktXgqbZyHojWOSHPxYOi8FEPW/42qJ9uttOfpMhCZDmbcRK+uCR9JJbVKmApdwfBItD1Fw1V5bW4tdu3bh8OHDUvrSQDM3K0KhEI4ePYoDBw5g3LhxcT+blDbeSY1kZGRo3HtSc9XV1VEziWPZV3w9Pr1xT8cSAcmb5CvG83a5J8YliV69cfVL5KRsUmqPg6tsru7M5o6ZNK2trfKayPbicTtOFiK/x+OBxWKRmanBYFAORenz7GnQO16ktMQKhzsqEZMBzzNJw+GO2potLS0yJKGXCHx4h69rw8fo+P8EIg0PbJKq0ROThwE4gfiqW/TgaTFKPt4HdI75kbQi6eL1euFwOOSC5RSdp0g5kYReCj7BNTMzEyNHjoTFYsGOHTukMZ+dnQ2bzYbW1la0trbKYiV6B+VESGli0XIcfGiCu8hNTU04duwYcnNzJfFigd5iMnSJMNwo5y68fuyOe18UftATg85Bwys8BEBhAzLq9eDDROFwWM6YIfuH1CnPuuCEpjFEcmToGKfTiZEjR0IIgaNHj8LhcGDw4MHw+/1obGzE0aNH0dbWhv79+0sSxouUJhYRgKadU3lpHgA8cuQIRo8erQknAJ22EJdaQKdnyMsNcWnHVRqPJemlGrfNAG0tVNr4ubljwM9F3hudiw80p6eny/O4XC6NQ0Ft07gf5VRZLB0rw7rdbjlyQcXfcnJykJWVBb/fj4yMDHg8HoRCIbn8y7d+XiE9JHpjqUAGoC0SGw6HpZdDy+XqR+i5UUuxJL66qNmsXaeQq1MKYvLZyUQSLh35FH29y04BUk5mUl00zEMeGnlyZBMRbDYbWlpaZN/dbrdUhVRnS1+TlQakeZZoMBiUWaUtLS1ynJKXqYx1DbFgEj0JTiQJvv76awwfPjzR3ThtUVVVhUGDBnV7TEpKrJycHADAoUOH4PF4Etyb5EVTUxOKiopQVVXVoxhUV6D4YDwJmSlJLDKoPR5Pr9ywbzsyMzN77T7F+yKndBxLIXmhiKVgCFKSWHa7HeXl5TGnhCl0IpH3KSW9QoXkR0pKLIXkhyKWgiFQxFIwBIpYCoYgJYm1YsUKucBQaWkpPv7440R3qc+QTKU6u0PKEeull17C/PnzUV5eju3bt2P8+PGYPn16VAWbbyuSqVRntxAphsmTJ4u5c+fKz+FwWBQWFoply5YlsFeJAQDx2muvyc+RSETk5+eLBx54QO5rbGwUdrtdrF27VgghxO7duwUA8cknn8hj3nrrLWEymcThw4d7rW8pJbGCwSC2bdumKUVpNptRVlYmS1GezjCqVOfJIKWIVVtbi3A43G0pytMZRpXqPBmkFLEUUgcpRaz+/fvDYrHIyZSEY8eOdVlm8nQCL9XJwe/PyZTqPBmkFLFsNhtKSko0pSgjkQjWrVsnS1GezujrUp3dotfcgD7Ciy++KOx2u1izZo3YvXu3mDNnjsjKyhLV1dWJ7lqfwOfziR07dogdO3YIAGL58uVix44dorKyUgghxH333SeysrLEG2+8IXbu3Ckuv/xyMXToUNHW1ibbuPjii8V3v/tdsXXrVvHBBx+IkSNHipkzZ/ZqP1OOWEII8dhjj4ni4mJhs9nE5MmTxZYtWxLdpT7Dhg0bBICobdasWUKIjpDDPffcI/Ly8oTdbhcXXnih2Lt3r6aNuro6MXPmTOF2u0VmZqaYPXu28Pl8vdpPlTajYAhSysZSSB0oYikYAkUsBUOgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMgSKWgiFQxFIwBIpYCoZAEUvBEChi9QBr1qyByWTCwYMHe/S7888/H2eeeWav9mXIkCG4/vrre7XN3oQiloIhUMRSMASKWAqGQBHrFPDGG2/gRz/6EQoLC2G32zF8+HD8z//8T8wVvABg27Zt+N73vgen04mhQ4fiySefjDomEAigvLwcI0aMgN1uR1FREe68884eLfSdDEjJlSmSBWvWrIHb7cb8+fPhdruxfv16LFq0CE1NTXjggQc0xzY0NOCHP/whZsyYgZkzZ+Lll1/GLbfcApvNhhtuuAFAx6zuyy67DB988AHmzJmDsWPH4t///jcefvhhfPXVV1FF1pIavTpL8VuO1atXCwCioqJCCCFEa2tr1DG/+MUvhMvlEn6/X+4777zzBADx0EMPyX2BQEBMmDBB5ObmimAwKIQQ4rnnnhNms1ls2rRJ0+aTTz4pAIgPP/xQ7hs8eLCcpJqMUKrwFOB0OuX/Pp8PtbW1OOecc9Da2oo9e/ZojrVarfjFL34hP9tsNvziF79ATU2NrKPwyiuvYOzYsRgzZgxqa2vldsEFFwAANmzY0AdX1TtQqvAUsGvXLtx9991Yv349mpqaNN95vV7N58LCQqSnp2v2jRo1CgBw8OBBnH322di3bx++/PJLDBgwIOb5UqkcpiLWSaKxsRHnnXceMjMzce+992L48OFwOBzYvn07fvvb3/ZoYW5CJBLBWWedheXLl8f8vqio6FS73WdQxDpJbNy4EXV1dXj11Vdx7rnnyv0VFRUxjz9y5Ihc+Jzw1VdfAeiIogPA8OHD8fnnn+PCCy+MWsU+1aBsrJMELQ8sWE2VYDCIJ554IubxoVAIq1at0hy7atUqDBgwACUlJQCAGTNm4PDhw3j66aejft/W1oaWlpbevARDoSTWSeJ73/sesrOzMWvWLPzqV7+CyWTCc8891+V6yYWFhfjDH/6AgwcPYtSoUXjppZfw2Wef4amnnpLrMv/sZz/Dyy+/jJtvvhkbNmzA97//fYTDYezZswcvv/wy3nnnHU1R2qRGot3SVII+3PDhhx+Ks88+WzidTlFYWCjuvPNO8c477wgAYsOGDfJ35513nhg3bpz49NNPxZQpU4TD4RCDBw8Wjz/+eNQ5gsGg+MMf/iDGjRsn7Ha7yM7OFiUlJWLJkiXC6/XK45I93KDqYykYAmVjKRgCRSwFQ6CIpWAIEkqs03kVr287Ekas030Vr287EuYVlpaWYtKkSXj88ccBdAxnFBUV4dZbb8Vdd92ViC4p9CISEiClVbwWLlwo9/VkFa9IJIIjR44gIyMj5Yc+UglCCPh8PhQWFsJs7l7ZJYRY3a3ipU83ATqyKnkG5eHDh3HGGWcY3k+F2KiqqsKgQYO6PSYlvMJly5bB4/HITZEqscjIyDjhMQkhVk9X8Vq4cCG8Xq/cqqqq+qqrCjEQj/mREGL1dBUvu92OzMxMzaaQ5EjUIOWprOLl9XpjLlSktr7Z+GB4V0hodsPJruKliJX8xErJ7IampiZ4PJ5Ed+O0hdfrPaE5khJeoULqQRFLwRAoYikYAkUsBUOgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMgSKWgiFQxFIwBIpYCoZAEUvBEChiKRgCRSwFQ6CIpWAIFLEUDIEiloIhUMRSMASKWAqGQBFLwRD0OrEWL14Mk8mk2caMGSO/9/v9mDt3Lvr16we3242rrroqqpyRQurDkIp+48aNw7/+9a/Ok1g7T3P77bfjH//4B1555RV4PB7MmzcPV155JT788EMjupISuOSSS+RCTQDw0EMPoa2tLYE96gWcWr2YaJSXl4vx48fH/K6xsVGkpaWJV155Re778ssvBQCxefPmuM+R6tVmzGazuPXWW8WmTZvEpk2bxNGjRzXXt2XLFrFp0yYxdepU4XK5RFpaWsL7zLeElDEqLy8XLpdLFBQUiKFDh4qf/vSnorKyUgghxLp16wQA0dDQoPlNcXGxWL58eZdt+v1+4fV65VZVVZXwm3uy28SJE8VvfvMbEQ6HT3gv29vbRTAYFM8884zIzMxMeN8TSqw333xTvPzyy+Lzzz8Xb7/9tpgyZYooLi4WTU1N4q9//auw2WxRv5k0aZK48847u2yzvLw84TfzRNvo0aPFQw89JJYsWdLlMdOmTROHDh06qfv69NNPC5PJlPDrTBix9GhoaBCZmZnimWeeOWliJbPEcrlcYu3ateLgwYNCCCFCoZB44IEHoo4bM2ZMXNUKu0I4HBZPPfWUcDgcCb/meIhleLghKysLo0aNwv79+5Gfn49gMIjGxkbNMV0VtSUkaw3SnJwcrF69Gtdeey0GDx4MoGPl1VhrN/v9/lPyfs1mM2666SasWLECbrf7pNvpKxhOrObmZhw4cAAFBQUoKSlBWlqapqjt3r17cejQoZhFbZMdw4YNw4wZM+I69uDBg1i7du0pn/OGG27o9iVMGpy0bO4Cd9xxh9i4caOoqKgQH374oSgrKxP9+/cXNTU1Qgghbr75ZlFcXCzWr18vVxydMmVKj86RLF6hzWYTDzzwQFT/XnzxxZjH33zzzaKtre2U7/HOnTvFgAEDkloV9jqxrrnmGlFQUCBsNpsYOHCguOaaa8T+/fvl921tbeKXv/ylyM7OFi6XS/zkJz+JcrdPhGQhFgBx9913R/WvK2IBEAsWLBChUOiU7zMtI3zaEKsvkMzE8nq9YvLkyV0e73A4hN/vP+V7UFlZKaZOnZq0xFJjhacAl8uFwsJCzb4vvviiT9ZdLC4ujukkJA1O+dVJAJJBYplMpphB3XHjxnX7O6vVKp599tleuQ9r164V6enpSSmxFLF6uNntdjFjxgyxcePGKFvpiSeeEC6X64RtTJs2rdfuRV5eXlISKyHLyqUazj77bDidTgDA888/j9zcXM3AOmHfvn1obW3ts3598cUXCAaDfXa+nkARqxuMGTMGs2fPxs9//nPk5OQkujtR+P3vf4+GhoZEdyMmFLFiYNasWZg1axYyMzMxdOhQAEB9fT3cbjdsNluXv1uwYAH+9re/4cCBA33V1aSFIlYMPPvss3juueei9s+dOxfnnXcegI41+y677DKNSiwsLMRFF12ElStX9llfkxa9ZkX2IZLFK1yyZElU39ra2sRNN93U5e+sVqt47bXXTvkefPTRR2LEiBEJuXblFRq82Ww2sWjRoqiAp8/nE1dffXXM3/RGgDQUCokHH3wwYdetAqQGIxgM4t5778Unn3yi2e92uw3NQDh69CgWLFhgWPu9AUUsg3DLLbegX79+UfvjWU/5RFiyZMkpt2E0FLEMwvPPPx8zFLBs2bJuPcvu4PP5cMMNN+DPf/7zqXbPcChi9QL+93//F5FIRLPP5/NF7QOAu+6666SCmoFAALfddhtWr14ds91kg1q6txdgNpsxbdo0jZrbvXs3jhw5EvPY2267DQ888AAsFktc7YfDYdx888145plneq3Pp4J4lu5VXmECNrPZLObNmyd8Pl+31xmJRMSmTZvE7Nmzk2YiBeL0ClWANAGIRCJ4/PHH0dbWhtzcXPzkJz/BpEmT5PfHjh3DH//4R0QiETz00EMIhUIJ7O3JQanCJEBxcTEGDBggP7e1tWH37t0J7FH3iEcVKmIp9BjxEEt5hQqGQBFLwRAoYikYAkUsBUOgiKVgCBSxFAyBIpaCIVDEUjAEPSbW+++/j0svvRSFhYUwmUx4/fXXNd8LIbBo0SIUFBTA6XSirKwM+/bt0xxTX1+P6667DpmZmcjKysKNN96I5ubmU7oQheRCj4nV0tKC8ePHY8WKFTG/v//++/Hoo4/iySefxNatW5Geno7p06fD7/fLY6677jrs2rUL7777Lv7+97/j/fffx5w5c07+KhSSD6eSZQBAMzEgEomI/Px8TWmfxsZGYbfbxdq1a4UQQuzevVsAEJ988ok85q233hImk0kcPnw4rvOmenZDqm99nvNeUVGB6upqlJWVyX0ejwelpaXYvHkzAGDz5s3IysrCxIkT5TFlZWUwm83YunVrzHYDgQCampo0m0Jyo1eJVV1dDQDIy8vT7M/Ly5PfVVdXIzc3V/O91WpFTk6OPEaPZcuWwePxyC2pq6woAEgRr3DhwoXwer1yq6qqSnSXFE6AXiUW1cbUF3HlxWvz8/NRU1Oj+T4UCqG+vr7L2prJWtxWoWv0KrGGDh2K/Px8TfHapqYmbN26VRavnTJlChobG7Ft2zZ5zPr16xGJRFBaWtqb3VFIJOJ0ACV8Pp/YsWOH2LFjhwAgli9fLnbs2CFXn7jvvvtEVlaWeOONN8TOnTvF5ZdfLoYOHaop6nrxxReL7373u2Lr1q3igw8+ECNHjhQzZ86Muw/KK0x+r7DHxNqwYUPMk82aNUsI0RFyuOeee0ReXp6w2+3iwgsvFHv37tW0UVdXJ2bOnCncbrfIzMwUs2fPPuHEAg5FrOQnlkpNVugxVGqyQsKgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMgSKWgiFISWKl4Lj5twrx3P+UJFZdXV2iu3Baw+fznfCYlKxBSku8HTp0SKXPdIOmpiYUFRWhqqqqV9K5hRDw+XxRyxXHQkoSy2zuELQej0flv8eB3pwnEO+LnJKqUCH5oYilYAhSklh2ux3l5eWw2+2J7kpSI5H3KSVz3hWSHykpsRSSH4pYCoZAEUvBEChiKRiClCTWihUrMGTIEDgcDpSWluLjjz9OdJf6DKlSqjPliPXSSy9h/vz5KC8vx/bt2zF+/HhMnz49qoLNtxUpU6oz7oIJSYLJkyeLuXPnys/hcFgUFhaKZcuWJbBXiQGQmFKd8SClJFYwGMS2bds0pSjNZjPKyspkKcrTGUaV6jwZpBSxamtrEQ6Huy1FeTrDqFKdJ4OUIpZC6iCliNW/f39YLJZuS1GezjCqVOfJIKWIZbPZUFJSoilFGYlEsG7dOlmK8nRGUpXq7DU3oI/w4osvCrvdLtasWSN2794t5syZI7KyskR1dXWiu9YnSIZSnfEg5YglhBCPPfaYKC4uFjabTUyePFls2bIl0V3qMyRDqc54oNJmFAxBStlYCqkDRSwFQ6CIpWAIFLEUDIEiloIhUMRSMASKWAqGQBFLwRAoYikYAkUsBUOgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIhlMJ544gmsWbOmT861e/duLF68GAcPHuyT83UHRSyD0dfEWrJkiSKWwrcYvTrn51uC7du3i4svvlhkZGSI9PR0ccEFF4jNmzfL78vLy0WsW7d69WoBQFRUVAghhBg8eHDUNK3zzjtPc+x7770n5syZI3JyckRGRob42c9+Jurr6zXtAhDl5eVR5xs8eLCc9kXt6bcNGzb0xi3pMVJyZQojsWvXLpxzzjnIzMzEnXfeibS0NKxatQrnn38+3nvvvR7NFn7kkUdw6623wu1247//+78BRBfsmDdvHrKysrB48WLs3bsXK1euRGVlJTZu3AiTyRT3uc4991z86le/wqOPPor/+q//wtixYwFA/u1zJITOSYwrrrhC2Gw2ceDAAbnvyJEjIiMjQ5x77rlCiPgllhBCjBs3TkqpWMeWlJSIYDAo999///0CgHjjjTfkPsQhsYQQ4pVXXkmolOJQNhZDOBzGP//5T1xxxRUYNmyY3F9QUICf/vSn+OCDD9DU1NSr55wzZw7S0tLk51tuuQVWqxVvvvlmr56nr6GIxXD8+HG0trZi9OjRUd+NHTsWkUgEVVVVvXrOkSNHaj673W4UFBQkhWd3KlDEOgl0ZfuEw+E+7Udfn68nUMRiGDBgAFwuF/bu3Rv13Z49e2A2m1FUVITs7GwAQGNjo+aYysrKqN+dyADXVzRubm7G0aNHMWTIELkvOzs76lzBYBBHjx7t0bn6EopYDBaLBRdddBHeeOMNjSo6duwYXnjhBUydOhWZmZkYPnw4gI7S2ISWlhb85S9/iWozPT09ihQcTz31FNrb2+XnlStXIhQK4ZJLLpH7hg8frjkX/U4vsdLT0wFEEz4RUOEGHZYuXYp3330XU6dOxS9/+UtYrVasWrUKgUAA999/PwDgoosuQnFxMW688Ub85je/gcViwZ///GcMGDAAhw4d0rRXUlKClStXYunSpRgxYgRyc3NxwQUXyO+DwSAuvPBCzJgxA3v37sUTTzyBqVOn4rLLLpPH/PznP8fNN9+Mq666Cj/4wQ/w+eef45133kH//v0155owYQIsFgv+8Ic/wOv1wm6344ILLoiqOdonSLRbmozYvn27mD59unC73cLlcolp06aJjz76SHPMtm3bRGlpqbDZbKK4uFgsX748Zrihurpa/OhHPxIZGRndBkizs7OF2+0W1113nairq9OcKxwOi9/+9reif//+wuVyienTp4v9+/dHhRuEEOLpp58Ww4YNExaLJaGhB1UfK0FYs2YNZs+ejU8++URTGvvbAmVjKRgCRSwFQ6CIpWAIEkqs03kVr+uvvx5CiG+lfQUkkFin+ype33YkzCssLS3FpEmT8PjjjwPoWAigqKgIt956K+66665EdEmhF5GQACmt4rVw4UK5ryereEUiERw5cgQZGRlJNYzxbYcQAj6fD4WFhTCbu1d2CSFWd6t47dmzJ+r4QCCAQCAgPx8+fBhnnHGG4f1UiI2qqioMGjSo22NSwitctmwZPB6P3BSpEouMjIwTHpMQYvV0Fa+FCxfC6/XKrbdzohR6hnjMj4QQq6ereNntdmRmZmo2hSRHQkYoxamt4uX1emPOSFFb32xer/eEzyih2Q0nu4qXIlbyEyslsxuamprg8XgS3Y3TFl6v94TmSEp4hQqpB0UsBUOgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMgSKWgiFQxFIwBIpYCoZAEUvBEChiKRgCRSwFQ6CIpWAIFLEUDIEiloIhUMRSMASKWAqGQBFLwRAoYikYAkUsBUPQ68RavHgxTCaTZhszZoz83u/3Y+7cuejXrx/cbjeuuuqqqHJGCqkPQyTWuHHjcPToUbl98MEH8rvbb78d//d//4dXXnkF7733Ho4cOYIrr7zSiG4YgpycHEyfPj3R3Uh6GFIq0mq1xiyg5vV68ac//QkvvPCCXKho9erVGDt2LLZs2YKzzz7biO70KlpbW6OWguspzGYzHA4HVq5cqVnJVY+FCxfi008/hd/vP6XzJQSnVogoGuXl5cLlcomCggIxdOhQ8dOf/lRUVlYKIYRYt26dACAaGho0v6FFjrqC3+8XXq9XblVVVQkv5dPVNnHixG6/N5vN4o477hDBYFBEIpFu72V7e7vYvn27uOiii0RmZmbCr422hNTHevPNN8XLL78sPv/8c/H222+LKVOmiOLiYtHU1CT++te/CpvNFvWbSZMmiTvvvLPLNmlx72TfLr30UvHll18Kk8kU8/tp06aJRx55RIRCoS6v9dlnnxUtLS1R+59++uku2z0tiKVHQ0ODyMzMFM8888xJEyuZJZbL5RJFRUVix44dora2VoRCIfHAAw/EPPauu+7SXNftt98uIpGI8Pl8IhgMCp/PJ6qqqmISLxwOi6eeeko4HI6EX3M8xDI83JCVlYVRo0Zh//79yM/PRzAYjFoBtKuitoRkrUGak5OD1atXo7KyEhMmTEC/fv1gsVhQVFQU1++XL18Ok8mEl156Cfv378fvfvc7CCFgsViijjWbzbjpppuwYsUKuN3u3r6UXofhxGpubsaBAwdQUFCAkpISpKWlaYra7t27F4cOHYpZ1DbZMWzYMMyYMeOUFzG48cYbMXbsWAwfPhxbtmzBN998g5dffjnmsTfccEO3L2HS4MTKrGe44447xMaNG0VFRYX48MMPRVlZmejfv7+oqakRQghx8803i+LiYrF+/Xrx6aefiilTpogpU6b06BzJUoPUZrOJBx54IKp/L774YszjMzIyolZq5Th69KioqakRlZWVorGxscvjdu7cKQYMGKBpOy8vTzidTpGdnS2WL18uCgsLv12q8JtvvsHMmTMxevRozJgxA/369cOWLVswYMAAAMDDDz+MH//4x7jqqqtw7rnnIj8/H6+++mpvd6NPEAwG4fV64z7e5/Ph1VdfjVoknLB//34cOnQIr776arc1Vs866yy5PjXhyiuvxODBg/G9730PCxcu1KwpnRDEJyOSC8kisQCIu+++O6pvkydP7vJ4h8MhWltbZaiBhxzo/8rKSvGXv/yl23tQWVkppk6dGtX+jBkzxKhRowQAw7zIpDDev81wuVwoLCzU7Pviiy+6XXfR7/fj2muvxZIlS/DRRx9h6dKl8rvnn38er7/+OgYOHIirr76623MXFxdHOQnp6en4+OOPceDAAQDAL3/5y7gdid6GItZJwmQyYenSpbjllls0++fMmXPC3/p8PrS1tcFsNsNms+Gzzz4DAIwZMwbDhw+HxWKB0+nsUX/MZjPmzp2L4cOH4+qrr4bZbMaKFSsStjxMQlb/SmWYTCbY7XYsWbIEv/rVrzTfrVy5EhUVFXG1c8MNN2Dx4sWora3F8ePHkZOTg0mTJp10v4QQePvtt7Fz50585zvfgUhw+X4lseKAzWbDtGnTMG3aNPznf/4nGhsbsWDBgqh40759+9Da2nrC9hoaGrB06VJccMEFePPNNzFo0CAsWrTolPoohMDOnTsBADt37kw4sZTE6gK33HILhgwZAqDDlpo3b16vtf3ZZ5+hoKAAf/zjH3HjjTfi17/+teb76upqfPrpp/jxj3/ca+fsc/TEG0sWGOkVjhs3Tqxfv15UVVWJurq6qK27gePDhw+L4cOHn/AckydPFseOHRM7duyI2U5bW5s4fPjwCe/DzJkzE+IJx+MVKomlw65du1BWVtbl96tWrUJ2djZMJhMuu+wyWK2dt7CwsBAXXXQRVq5c2e050tPTkZubiyeffBLjx4+Pitw7HI4ob1OPysrK5F638YTUS0IkQxzLZDKJ66+/XrzxxhuavrW1tYmbbrqp299Omzaty2tbtWqVWLJkyQnvwerVqxN27UmR3WAEkoFYtHk8HvHJJ59o+ufz+cTVV18d83i73S7z044dOyY2btwo/va3v4lQKCT+8pe/iJdeeumEajAUCokHH3wwqYmlvMJThNfrjcrwdLvdXWYgmEwmucj6c889h3PPPRfZ2dn48ssvsWnTJsyYMeOEavBPf/oTFixY0DsXYBCUjZVA3HHHHQAAp9MZM1WmK+jHCZMRSmIlAUpKSpCVlYVBgwZ1e5zP58MNN9wQdxA2kVDE6gW88MILePDBB/Hpp5+edBsFBQUoLy/v8vtAIIDbbrsNq1evRiQSOenz9BlO3oROHJLJeOfbiBEjRFlZmSgrK+syH8psNovbb79dph9v2rRJvPnmm91ebygUEj//+c8Tfn20Ka8wSTez2SzmzZsnPvroIxEKhbqcXBGJRMSmTZvE7Nmzk2YiRbzEUouNJxAOhwO33347AODOO+9EVlYWAOD3v/89mpubEYlE8NBDDyEUCiWwl9GIZ7FxRawkwZlnngm73Q6gYyyxqyzTZEA8xFLhhiTBF198kegu9CqUV6hgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVD0GNivf/++7j00ktRWFgIk8mE119/XfO9EAKLFi1CQUEBnE4nysrKokor1tfX47rrrkNmZiaysrJw4403orm5+ZQuRCG50GNitbS0YPz48VixYkXM7++//348+uijePLJJ7F161akp6dj+vTpmizL6667Drt27cK7776Lv//973j//ffjmkGskEI4lSwDAOK1117TjMbn5+drSvs0NjYKu90u1q5dK4QQYvfu3QKAJk/8rbfeEiaTKa4pT0KkfnZDqm99nvNeUVGB6upqzfQpj8eD0tJSbN68GQCwefNmZGVlYeLEifKYsrIymM1mbN26NWa7gUAATU1Nmk0hudGrxKqurgYAOVmAkJeXJ7+rrq5Gbm6u5nur1YqcnBx5jB7Lli2Dx+ORW6IqqCjEj5TwChcuXAiv1yu3pJ6oqQCgl4lFtTH1S5jw4rX5+fmoqanRfB8KhVBfX99lbc1kLW6r0DV6lVhDhw5Ffn6+pnhtU1MTtm7dKovXTpkyBY2Njdi2bZs8Zv369YhEIigtLe3N7igkEnE6gBI+n0/s2LFD7NixQwAQy5cvFzt27JCze++77z6RlZUl3njjDbFz505x+eWXi6FDh4q2tjbZxsUXXyy++93viq1bt4oPPvhAjBw5UsycOTPuPiivMPm9wh4Ta8OGDTFPNmvWLCFER8jhnnvuEXl5ecJut4sLL7xQ7N27V9NGXV2dmDlzpnC73SIzM1PMnj1b+Hy+uPugiJX8xFI57wo9Rjw57ynhFSqkHhSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMgSKWgiFQxFIwBClJrBQcN/9WIZ77n5LEqqurS3QXTmv4fL4THpOSFf1ycnIAAIcOHVLpM92gqakJRUVFqKqq6pV0biEEfD7fCVfOAFKUWGZzh6D1eDwq/z0O9OY8gXhf5JRUhQrJD0UsBUOQksSy2+0oLy+X5asVYiOR9yklc94Vkh8pKbEUkh+KWAqGQBFLwRAoYikYgpQk1ooVKzBkyBA4HA6Ulpbi448/TnSX+gypUqoz5Yj10ksvYf78+SgvL8f27dsxfvx4TJ8+PaqCzbcVKVOqM+6CCUmCyZMni7lz58rP4XBYFBYWimXLliWwV4kBkJhSnfEgpSRWMBjEtm3bNKUozWYzysrKZCnK0xlGleo8GaQUsWpraxEOh7stRXk6w6hSnSeDlCKWQuogpYjVv39/WCyWbktRns4wqlTnySCliGWz2VBSUqIpRRmJRLBu3TpZivJ0RlKV6uw1N6CP8OKLLwq73S7WrFkjdu/eLebMmSOysrJEdXV1orvWJ0iGUp3xIOWIJYQQjz32mCguLhY2m01MnjxZbNmyJdFd6jMkQ6nOeKDSZhQMQUrZWAqpA0UsBUOgiKVgCBSxFAyBIpaCIVDEUjAEilgKhkARS8EQKGIpGAJFLAVDoIilYAgUsRQMwf8DprZY9FtCpxEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best IoU: 0.5263, Best Dice: 0.6666, Best Acc: 0.9622\n"
     ]
    }
   ],
   "source": [
    "# -----------------------------\n",
    "# Threshold method application\n",
    "# -----------------------------\n",
    "\n",
    "# empty the cache for model and training process\n",
    "torch.cuda.empty_cache()\n",
    "\n",
    "# Do the segmentation to exteract the spine part\n",
    "# Threshold segmentation hyper parameters\n",
    "minThreshold = 250\n",
    "maxThreshold = 1500\n",
    "\n",
    "# the variables for evaluation\n",
    "total_correct_pixels = 0\n",
    "total_pixels = 0\n",
    "total_iou = 0.0\n",
    "total_dice = 0.0\n",
    "\n",
    "# set the hyper parameters and the paths 设置超参数和路径\n",
    "data_path = \"./dataset/Spineweb_dataset15/test_1/\"    # the path of the dataset\n",
    "file_list = os.listdir(data_path)  # list all files and folders in the folder(data path)\n",
    "file_count = len(file_list)  # get the numer of files\n",
    "print(\"Length of the dataset:\",file_count)\n",
    "resize_scale = (128,128)    # set the resize scale for images\n",
    "print(\"----------------  Loading the dataset and evaluate  ----------------\")\n",
    "data_count = 0\n",
    "poster_threshold = 200\n",
    "for index in range(file_count):     # iterate the whole dataset\n",
    "    metadata = np.load(data_path+str(index)+\".npy\")     # load the data\n",
    "    image = np.array(Image.fromarray(metadata[0]).resize(resize_scale))     # Using the PIL to resize the transfer into npy again\n",
    "    label = Image.fromarray((metadata[1] * 255).astype('uint8')).resize(resize_scale)     # Using the PIL to resize the transfer into npy again\n",
    "    label = np.array(label)\n",
    "    label[label < poster_threshold] = 0\n",
    "    label[label >= poster_threshold] = 1\n",
    "    \n",
    "    # print(type(image),image.shape,image)\n",
    "    # print(type(label),label.shape,label)\n",
    "    output = CTimageLIB.thresholdSegmentation(image,minThreshold,maxThreshold)\n",
    "    # evaluations\n",
    "    current_correct_pixels,current_pixels = calculate_pixel_accuracy(output,label)\n",
    "    current_iou = calculate_iou(output,label)\n",
    "    current_dice = calculate_dice(output,label)\n",
    "    # add the results to total one\n",
    "    total_correct_pixels += current_correct_pixels\n",
    "    total_pixels += current_pixels\n",
    "    total_iou += current_iou\n",
    "    total_dice += current_dice\n",
    "    data_count += 1\n",
    "    # visualize the results\n",
    "    if index == 100:\n",
    "        visualize_results(image,label,output)\n",
    "\n",
    "final_acc = total_correct_pixels/total_pixels\n",
    "final_iou = total_iou/data_count\n",
    "final_dice = total_dice/data_count\n",
    "\n",
    "print(f\"Best IoU: {final_iou:.4f}, Best Dice: {final_dice:.4f}, Best Acc: {final_acc:.4f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Length of the dataset: 3394\n",
      "----------------  Loading the dataset and evaluate  ----------------\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJYAAAHWCAYAAACVLnA6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABTnUlEQVR4nO29e3iU1bU//plLMpPbzCQQEqKBcFW8YR+QgMWKiqJVrMpzVMrxQcsRr/S01KqcbxWxtnhpxaqo2PZIryrYqk+reIuoVQEVovagIiiXCIZLrpNMMpPM7N8f+X121ju5MIG8mRnZn+eZJ5l33st+3/1511p7rbXXdiilFAwM+hnOZDfA4JsJQywDW2CIZWALDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkOsbrBy5Uo4HA7s2LEj2U1JWxhiGdgCh4kVdkU0GkVbWxs8Hg8cDkeym5OWMMQysAVGFXaDeBurrKwMF1xwAd544w1MnDgRWVlZOPHEE/HGG28AAP7+97/jxBNPhNfrxYQJE1BZWWk538cff4wrr7wSI0eOhNfrRXFxMX7wgx+gpqamy7V5Da/Xi1GjRmHFihW44447upWcf/7znzFhwgRkZWWhoKAAl19+Oaqqqvr9eRwSlEEXPPHEEwqA2r59u1JKqeHDh6tjjjlGDR06VN1xxx1q2bJl6qijjlK5ubnqz3/+sxo2bJi6++671d133638fr8aPXq0ikaj+ny/+tWv1GmnnabuvPNO9fjjj6v//u//VllZWWrSpEkqFovp/TZt2qQ8Ho8qKytTd999t/rFL36hSkpK1Pjx41V8V911113K4XCoyy67TD3yyCNqyZIlavDgwaqsrEzV1dUNxGPqFYZY3aA7YgFQ7777rt7n5ZdfVgBUVlaW2rlzp96+YsUKBUCtXbtWbwuFQl2u8eSTTyoA6q233tLbZs6cqbKzs9Xu3bv1tq1btyq3220h1o4dO5TL5VK/+MUvLOf897//rdxud5ftyYBRhQniuOOOw5QpU/T38vJyAMCZZ56JYcOGddn+5Zdf6m1ZWVn6/9bWVhw4cACTJ08GAGzatAlAx4Dhtddew0UXXYSSkhK9/+jRo3HeeedZ2vL3v/8dsVgMl156KQ4cOKA/xcXFGDNmDNauXdtft33IcCe7AekCSR4A8Pv9AIDS0tJut9fV1elttbW1WLJkCZ566ins27fPsn9DQwMAYN++fWhpacHo0aO7XDt+29atW6GUwpgxY7pta0ZGRiK3ZCsMsRKEy+Xq03YlBtuXXnop3n33Xfz0pz/FySefjNzcXMRiMZx77rmIxWJ9bkssFoPD4cCaNWu6vX5ubm6fz9nfMMSyGXV1daioqMCSJUtw++236+1bt2617DdkyBB4vV5s27atyznit40aNQpKKYwYMQJjx461p+GHCWNj2QxKFBXnLnzggQe67Dd9+nQ899xz2LNnj96+bds2rFmzxrLvJZdcApfLhSVLlnQ5r1KqWzfGQMNILJvh8/nwne98B/feey/a2tpw1FFH4ZVXXsH27du77HvHHXfglVdewbe//W1cd911iEajePjhh3HCCSfgww8/1PuNGjUKd911FxYtWoQdO3bgoosuQl5eHrZv345nn30W8+fPx0033TSAd9kVhlgDgL/+9a9YsGABli9fDqUUzjnnHKxZs8Yy+gOACRMmYM2aNbjppptw2223obS0FHfeeSc+/fRTfPbZZ5Z9b731VowdOxbLli3DkiVLAHQMJM455xxceOGFA3ZvPcGEdNIAF110ETZv3tzFLktlGBsrxdDS0mL5vnXrVrz44ouYNm1achp0iDASK8UwdOhQHVfcuXMnHn30UYTDYVRWVvbot0pFGBsrxXDuuefiySefRHV1NTweD6ZMmYJf/vKXaUUqAMkNQj/88MNq+PDhyuPxqEmTJqkNGzYkszkG/Yik2VhPP/00Fi5ciMWLF2PTpk0YP348ZsyY0SXkYZCeSJqNVV5ejlNOOQUPP/wwgI4wRWlpKRYsWIBbb701GU0y6EckxcaKRCLYuHEjFi1apLc5nU5Mnz4d69atO+jxsVgMe/bsQV5enkkdHkAopRAMBlFSUgKns3dllxRiHThwANFoFEVFRZbtRUVFXRyBABAOhxEOh/X33bt347jjjrO9nQbdo6qqCkcffXSv+6TFqHDp0qXauywxduxYtLW16bdHKQWn04lYLAan04n29naLRHM4HHA4HIhGo5bvTqcTDocDSilEo1G43W54vV7k5OQAANrb25GZmYlYLIZYLIbMzEy0t7fD6XQiGo3C4XCgvb0d0WhUxwZbW1vh8Xjg8/ng9XoBdEhat9sNh8OBtrY2KKX0OaPRKDIyMuD1ehGNRnUGQ3Nzs96P9+n1epGRkQGlFNrb21FTU4OmpiZ4PB74/X59z5mZmfq67e3tcLvd+lyRSAQulwuRSEQ/O4/Ho6/jcrmQlZWln6HD4UAsFsPrr7+OvLy8g/ZZUog1ePBguFwu7N2717J97969KC4u7rL/okWLsHDhQv29sbERpaWlcDgcuiNdLhfa29s1udzuzltTSukHw47l/yQS4XK5kJGRgezsbLjdbsRiMU0Mh8MBt9utH7wkqsfjQTgchsfjgcvlQmZmJtxuNzIzM5GRkaEJ73A4kJGRocmdkZGhO55k9Xq9aGlpQTgc1mTiMTze5XLp9ns8HkQiETidTmRmZkIppdsXjUbR3t4OAPo3kowpO5FIBG63W7ezra0NWVlZ+ny897a2Nv3/wZCUUWFmZiYmTJiAiooKvS0Wi6GiosKSpUnwzZcfHkNJwocppYDsTHZQNBqF0+nUH+4PdNh5Ho9Hk0+SEICFuJmZmVrSyfuitOH5eQ55HUoJAAiFQmhubkZDQwNaWloQi8XQ1tamJSDHVlIq8574cvC67HySn9f3eDyamHxOJBEAy/eMjAxkZWUhHA5rEobDYd2eRJE0Vbhw4ULMnTsXEydOxKRJk/DAAw+gubkZV111VZ/PRQJRMvGh8+ECsKgS7stj2fFAh5rh7+3t7cjIyNDShNJNdq7L5YLT6UQkEoFSSqsotiFesjkcDkQiET13MRKJwOFwwOv1IhwOo7m5Wd8XVTPJSyLx/FS9Xq9XS5yMjAy0tbVpsrndbrhcLn1PbJN8CeNfPIfDoaWlUgput1tfN1EkjViXXXYZ9u/fj9tvvx3V1dU4+eST8dJLL3Ux6A8G+RZSUrEjehPZ3N/tdmsbiuKfH6oQt9utVQ2vScnAjgOgiSizOiORiLZV2JkulwttbW2W7FFKCw5SpKr2eDyWdmdmZupzSeJxP74M8iWQLxq38b4zMzMRiUQQi8XQ0tKipVw4HEZrayt8Pl+PmbI9IanG+4033ogbb7zxkI+Xb5AkieqYfWRJsosnmTRS2UkZGRmaKPGSgqpEEhKAtkeADmLxO0lD4rNDgY6Op9SS55QE4D1wgMB20GYiYSmVqHLZDt4TyUQ1SRIppbRKdjqdCIfDaGpqAgAEAgF4PB79THgffUHaZzdI0rDj2FHsNGmrxB9HAsmBAB+ifEvlSJJ/44kt1R7Vc1tbmz6W7WDHZ2VlaalEskYiEbS2tur9KTlJuLa2NoTDYUSjUYTDYW1Peb1efU0A+iXjcyAomVpbWxGLxRAKhdDQ0IBQKKRfDP5G9c9BRV+Q1sRiJ8ublqMd7iNHifI3qRKkLcIOknYZj6V6pCRgJ8lOlW2jdCEJpQukOxLwmiRhKBSyjCZp/7D9PBdtPanqeM9y4CBHtW63G36/Hzk5OXA4HHoE29bWhvr6ejQ1NWmp1ldVmNbEih/RSfVBVcDfpErisZQucnjODztbjr7kaEyqUr7pHI4DnSTkteOlilJKE0QSPxAIIC8vD/n5+cjLy9O/ZWdnW/aVtiRVHd0JcjQsCS5NBT6PtrY2+Hw+DB06FD6fD5mZmVptUpJmZGRoKZoo0sJBejBIw5SdLiWEVFlutxttbW0WUkknKf/nMdK+oU+Hv5OA3E+ORkl4SjS3260lgxzyUzUqpbTEAKD/p93ndrsRCoW0BCFxSBLagFK1Uj3zvjnCo/SLxWIIh8OoqanB4MGDUVhYiH379qG5uVnbf9Kg74vUSmuJJW0rdiBtJXYG95NvL3+XLgSqQUofAJaRIb3r0lAm6Xgd6feSpHM6nWhqarLYKpJgssNIJto2bGcoFNL+PLYf6Ah3kWxsEyUpt3GkKrfJZxQOh9HY2AiXy4VBgwbB4/Hol8fpdKK+vh6tra2WsNrBkNbEkh0jXQTSWKVKlFJLGt48jsTi/iSB9N/QeKbxzYdPcnL/lpYWiy3GEEooFLL4zSiNvF4vvF6vVkNUn3RJ0E8ViUTQ3t6OrKwsLdHoumCIhvfHl4Vtly8CyU1Vn5GRgVAohKamJni9XmRnZ3fx7zESkCjSmlhA58iJKkiOyKSqkJAGPf+XXmjaGJI4ctTI60pVw06SKoPqhu1paWnRI1WHw4Hc3Fzt8c7OztaF3qSU4chSSkMpkQDoUSLPTcJK4573w/igJAm9901NTRb/HNtDd0RjY2PC/ZLWxJLGOCUG0DnairefuI2edGknxXvo6aaIvx7fYKoiSpHW1laL64LnkjZTKBRCMBi0SEbpX5O2GSWiVHtUk9ID39zcjFAoBKCD7FSLbJdUuSRpOBzWkp6VC+lbo7SVQXU+m75IrLQ23uUb3J2fih3IB9OdRAM6/V5AJ7GkHSXJSBJK24vH8cNrUdrJ7XRC5uTkaA86O5jtyMnJsXSiDA+5XC60trYiFAp1CQHxXmSwmP+73W7tn4pEIsjOzraMYqlO6+vr9f5Uz7FYDI2NjV1mEPWGtCYW0GlgS9XAbdKuiD9GqjGZ4SDdASSRHH63trYiNzdXB6wZU8vPz7eMDAFYSM/RaHt7Ow4cOID29nbk5eXpbAiCMT7aTPJ8bW1taGlpQVtbG2pqaixxRtlWHi9Hf0zT4QtB90F7e7tOw2lvb0dLS4s+36BBg5Cbm4twOAyv14v9+/cn3C9pTSypdqR0kUFVoNPIBzo6RwaK4x2XdCnQZgFgkXr0jOfk5KClpcVybo/Ho10ZPD/jcfRjNTQ0oL29HdXV1WhsbERWVlaXTAOqI44IafOw45lxAHS6JZhLxeuSlMwd43OgTUpQpdIGi0ajaGlp0b4rjpRbW1stWRkHQ1oTSz4kvtUkGJ2ITqdTuyG4TdocACx2Gu0KSbx4VwUlB4+RoSEGkuO96Tw2OzsbkUgEbW1taGxs1J0XnxMmXScAtKHvcrk0qaRnXap3+tyky0DacvK+6M1n6Ic5ZZmZmWhra0NTU5N2ffQluyGtjXcAFgOTtpGMr0nHpgwM87sM0soRoFQtkngej0dnK3A0JcNFlA5Ap00nk/noKM3Ly9Nk4bVo84TDYU0E5kU1Nzdrqcr6VzJITVXXXWyUajxebXIE29raqglEyVpaWor8/HyLv6u3bJF4pLXEIuQQnqqIkE5JKaHkPkopnYVJqSXtMBJRGvD09chRIkeNXq/X4k+T+0mS+v1+HdohAd1utzbIY7EYPB6PlhhNTU3a10WJSGJ4vV4dXpL+LUpfZjawLVRv/PDeAehEyv379yMUCmkpmRb5WP2B+DeIdhGT9eI7oKcHQwJwBARYDW8AFgJRClIiyhEiO02qZjl8JwFkdgXTbkiiQCCgz0kCMM2FNhvzw6jC6urqUFhYCKDTkStfFD4DjlIbGhq0JJN5Z/S3tbS06GgB778vGQ5pTyxKH0oSABa7RsbyAKunPT7thR1FY1oGYeXQXHrHadBSXfDclB4kC2DNOZeOWXltr9erBwdtbW16Ukc0GtWpLTJsQynJfKrs7GyLOpSOW8CawyalKEkVi8UQDAa14U83hcvl6lMJyrS2saQjlPaRNMrl8FpmUwKw2EIyYEspEp/DJPfn8XR+hkIhhEIhtLS0aMLJgYUMcNNWkSk73If7U93ReGfWaHZ2trbLAHRpa0tLi8634vX50lHNUsozHihJxWcon5l0l/SlaG5aSyygM3NUppJIm0qOzuirkXlT8TNYpDNSxtmkFJKpKF6vF0opNDU1aVLFq0QZt+P1pO0m3QRSbZM4MhOW6dNAhy+KqpESidPAaLPJZyHzwmTiYjQaRU5OjrYzqbqbmpr0s6LRnyjSnlhAp8SR+VjsQD6YSCSinZHxoZTu7Aip7uJ9XtyXJMjIyIDf79fpK7I9ra2tOlNUhmJILBkuoYSgO8Hj8ehO5nXa29st8x0B6PmINLAbGhqQk5OjjXYAOotB2lkkisfjQXFxMQKBAJqbm9HW1oa9e/fql4XRgL6MCtNaFcq3UXrRJUHkKEzaHJRM0msv02BkMhxVirwOpRjPn52drTs+HtKLz2PZ6VJ6SruQRjdfFpnJKfPcZVyTqo7uA04n48iRWRHS3SDbHgqFNEkpoThQoZmQcN8kvGcKIt59IAPNfKDShgGs3nq5Pd6Q518ZGJb/S7DjSEgZa2S6i3RZOJ1OnbTHADJdAbRpZJinp/byOpRqJARJEK/K6evyeDxa6gUCARQWFqK1tRVff/01GhsbUVtbq10YTmfnJNgjxt1AtSYDytJA7y5dhvtIfxW/M9epO+Nbkqu5udkya1jmhHFyp0w0ZN47jV8GkaWRLKfNc5QniSZHvNJxSwLR3pJrLFId838+J5liM2jQIPh8Pj3w4Ox0vpQkFm3HRJHWxAI6O10ao1K1dacuAauEYgfFh3lk/JEdweE+PeAMvfB3Elbmc8lANn1L7HSSjd+lcc82sQ1yX/q3GA1gPhWP5QQPAFr98T7YBim9/X6/jjkyc4KjSU5n607N94S0JxYhsxLkSIpvONCZxSBDN0CnrwuAxf4iWSVBpT+H0srhcOjgNFVRLNY5e4ehlLy8PNTX1+sgL6UjSZOfn49AIKAlVjAYRGtrq8UZ63R2zLpm9oOcUk+PPwlK45/PRqbeUH0Gg0F9jqysLB1IZ1IiY4lA10UQekNaE4uinSqRnRwfMJW2knSkkiycnCnjhZKYlEQyXZnXkFkPUnVIRym92bymtH8yMzNx1FFHYcyYMRavPEd3rOvAnHN66qWBT1UrXSZUX1TZHLTE24gHDhzQzzA/P19L0+7CUUeMKuRNS5UlbZGeIN0NWVlZmkQyhiilVXyaDc/NTAZKEA7pZdxNZnzKCRKtra3Iy8tDWVkZysrKdHIf86hk+AfotCdbW1t1LJB2W0tLiyV7g7lWbKvcjy8D7UOOit1uNw4cOKBfLN4/sy/48iXcNwnvmSDu+P+XmZWfY489Vv/e2tqKG264QSeRzZo1q0s5o75AhiVog3TnyKOalCotPn8J6IztkTysPiP9WdyHHROLxZCXl4eCggI9F1BKNOlqYH6Ux+PB4MGDccwxx+iaVuxUGvdUwXLChSSarGhIqSVfqHiXBSWZrInl9XpRWFiIESNGoLS0FF6vV//Oe2NyYXdJkz3BFol1/PHH47XXXuu8iOjoH//4x3jhhRewevVq+P1+3HjjjbjkkkvwzjvvHNK1pHiXKk9mJwBW5yfJyM5gsTMZaKbhKgcFlIjdee7Hjh2LSy+9FB6PB2vXrsWbb76p7RpKKKCTqJRWLpdLz94hotGoVqfSe856DXRf0PZyOp16+pacWEqJ1dTUhLy8PH2urKwsi49ryJAh8Pl8WvIGg0Ftw/F4qvZEYQux3G53twXUGhoa8Pvf/x5//etfceaZZwIAnnjiCYwbNw7r16/Xq472Bd1lLUiDXSbo8TdJEqAz7YYdReOWEomuAKbkxA/5gY4VWCdPngyPxwOv14svv/xSr8LFN166IeiU5GQIn8+nVaDT6URLS4tlMMDgtEwAZNDa4XAgOztbjxIBq7qX2Q1UzQQHFkCn1JfxQgB6EML9EoEtDtKtW7eipKQEI0eOxJw5c7Br1y4AwMaNG9HW1obp06frfY899lgMGzas16K2nHokP0Bn8Fm6FOJJxtFPfEaBjOeRLNJ1wIfMgDDQdeVSuV9OTg4aGxtRX1+PkSNHYurUqcjKyrIk2dFZCkBnL1BasY1tbW1oaGhAY2MjmpubsW/fPp3GTMemzI7lPefm5sLv9yMzMxM5OTkYNWoUTjjhBIwePVqn09B1QZ8VJen+/ft1KKe2thYAtLdd+uLkqrEHQ78Tq7y8HCtXrsRLL72ERx99FNu3b8dpp52GYDCI6upqZGZmIhAIWI4pKipCdXV1j+dcunQp/H6//nC5XEbsAVgII9WKdHLKoK+0q+iR5oe2hRwdyn3k9TweD4qKilBYWIjq6mq9XNzUqVO1DZSZmQmv16szCoDO3DGSif9LcktjPBgMavVM94SUzE6nE3l5eSgqKkJZWRkKCwvh8/lQVFSE0tJSFBQUQCmlp8tL7NmzB1999ZWWsBxd84WSI+pE0e+qUC6MfdJJJ6G8vBzDhw/HqlWrLItu9wU91SCl/SJtLAAWVQFYSzPSdiEoLZilICdQyECxnHZPYrrdbvh8Pnz729+Gx+PB1q1b4Xa7NfGl70mmSnNun1StsmYWAD2ZgRKZGRS8rizQRukXDAYxePBg5OXlWWb4eDweBAIBhMNhtLS0wOPx6CJtlJT79+9He3u7tsVqamq0jcbrJ914lwgEAhg7diy2bduGs88+W89dk1Krp6K2hHwQEvEed+l6ADodonJURiOUwVqmtEgpIW0s6ZSkxMjMzMTw4cNRVlaGo446CgUFBdixYwccDgeKiooQDofxwgsv6Bk5MvwSn/vEUaIkDUsXSRVOSSYlMsnIewU6J3dIXxlfQLoMmECYl5dnsSVlmvPQoUNRV1eHYDBoKbKbKGwnVlNTE7744gtcccUVmDBhAjIyMlBRUYFZs2YBALZs2YJdu3Z1W9T2YJDJeUCn15yjOsA6g4VBVTkRglIE6By9MkuB/iin04ni4mLk5+ejuLgYRx99NEaPHo22tjZUVVVh165d2mk5bNgwfPbZZ6isrNT+JZkmw2vRa09JSF8TJSRzpygtuT0ajVpmLtOlED8C5nepLmWSIUeYNMilSREKhZCbm4v8/HzLy8tJKYmg34l10003YebMmRg+fDj27NmDxYsXw+VyYfbs2fD7/Zg3bx4WLlyIgoIC+Hw+LFiwAFOmTDmkEaH0CvPhshM5WpRuBBKKxrQMHlMdMa7GWlFHHXUUSkpKUFBQoEdmDQ0N2LVrl8U/lJeXh7y8PGRmZuosAUoHEpjkkBkMDDFRYtKhSVLK9BbafNKVIj3+tMcYOwQ6E/mam5u1gU4Hq5R8JBttP6Y5+/1+hEIhS+pOIuh3Yn311VeYPXs2ampqUFhYiKlTp2L9+vV6ZLJs2TI4nU7MmjUL4XAYM2bMwCOPPHJI15IuBUkkqjR2BAPAMsgr87CoJo499ljk5OTA5XLpst8ejwdNTU3YvXu3LpohZymzaP/gwYPh8/mQk5Oj6zMAnTYdj5Geb3ZuvK0V7yogiXgMScr9qGaZXx8MBnWdeqDDXbBv3z40NDRY6m7JnH7mXbGddJEwh6uxsbFP6rDfifXUU0/1+rvX68Xy5cuxfPnyw75WvPhnp/C7dIhKj7mUAEx1ycnJwbBhw+B2u7Fv3z49EqPB29DQAKUU8vPzMXjwYL2eTFNTkyXBTgZs5SgUgEUNAdC2Dn+XZM3OztbSi+qaFY4ZZpHuFkqiSCSC3bt3Izc3VweU29raUFdXB4fDgUAgYJGabJ8cEXPihlShcjCUCNI6VijVHEGbgg9Oxva6C+mwQ2KxGCorK3X6CbMwAevEi6+//lqnIp9wwgkYOXKk7hhKnhNPPBHvvfce9u7dayGY9Ow7HA5ty9FxSXICndVp2KmyVBHQGXGQkYfm5mbU1dWhvb2z6L9Mt5FmgpSE8lwOh0PXHqWfK94HlwjSmlgy2i7DDfFvlvTDUHXER+yj0Sh2796tg7+UbjLNhHE8km/v3r348ssvcdJJJ2HMmDHIysqCUgrjxo3DtGnT8Nxzz1kmX5D0JCBtJGnfORwd6+cww7O+vl4TgSosIyNDS1Deh8vl0pMfcnJy9OhOqn+WgIw36qXpIFOSY7GOUpK85hGTjyWndHWX8Sn/Ssem9CZTuvEcnPXC7dIIpqqRtkZVVRVqa2vhdrsxceJEDBo0CNnZ2Tj//PPx8ccfa98WR1tsN8nKUBJHZ1RfMrZIVc1RLSWpdClINwLDRZSMrMHg8Xi0cc/fpBrkC5Cbm6sry0ib9YhRhdKolRkJlAYU7fItlT4knoNvpnRdSIcmSRuLxbSTV8bTgsEg3n33XYwaNUoThfFAqjEm09EgpppiDSr5UlB90QnKAQX9aqFQyGLLUY1lZ2frMBKvzUkUHo9He9bjXzYa8jI7xOv1orm5GS6XS7s3kh6EHijInHUJaVeReDKeKIOyDLnwfzmaBKDVlsz0lHYcCSNjfQDw73//G3v27NFStampSTsaHQ6Htpk4epMpy06nEzk5OZYVL2j3NTc3W5zAbCPQOfdPqizmwCulkJ2d3W0SI5+hlEh8QejEpeRMFGk9S0fG/+JTYgg+MEqp+DgZVQAzGGT4giqDpJOOS5lvFgwGAUAfv3//frzzzjuor6/XUobnknYNjXO/3w+fz6djjwD0wKK5uRnBYFDnoXO0Ju+XhJfBcrpXpGHOF4jSmAQlMfks6QeTvjWeL1GktcSSdkF88JkPkPvRKSr9VzLZTzofZQYozy3/53mYbhyLxVBXV4eNGzfC4/Hgiy++QGVlpS6sRuNcVoEJh8PYt28fhg4dipycHB0bZJ5Uc3OzXrewtbVVZ0LQCcrz0N1BScp7kjXl+XLxOHlvMu2Iz0BOCqHLg5MwEkXaEyvecy1ntci3Feh8iNKZKFNogM64GtC5+ICsr0V7hB92Qm1tLV577TWsX78ewWBQqz7aUTJ4zZEaR5Vut1svDiojAjTiSSy+KCQ2Dfz4zNl4aUPVy2fGNsggPs8lncaAtQR5SgWh7YR8G6XRTX+VdIayM+V3aVOQYHIYTmdmfBYF/UGyAFskEkFTUxPq6+v1xIZBgwb1+qa3trZi165daG1tRUlJCQoLCy0DjAMHDmgPOKfgE5xgwXQZOZFCzprmNqp0WfGGpOF2GaFgSUpZ2umIcTfEp/NKyOCzzCiQPiXAuuSvPJbnlwXNpHsCsJYAlxWJac9wRg1fAI/Hoyvz8ZwtLS3YtWsX9u3bhyFDhiAQCGg7jvFGOQJVSul1BKmypSQj6SnV6H+S0l2OAKW0Zrvb29vR3Nys3S3x0iwRpD2xKJlk3DDe7oofIZJkhHy4fEPZadJPlJ2drSvcyURBp9Opl7vlOTgpQdYq5UiQbeG8P15/3759aGxs1IQArKlBfDHoTOVcxpaWFu3nam5u1qpbumLkS0E1Ll0cbCMAXSOVsUkG36WtedC+SXjPFIQc0cghNADLW0zEpynLTpNDbukQjA8cUxJRPfF6dEjm5uYiKyvLopqlK4JSkBNWuY/8n23gaAyARXJxJEmnKlOieX65tg6n3jNDQRazpcSlRKUtxcwKDiIo7ZKaNpMMdOd2kLaUdA9wdEhQlTBmKGfm0O5ih/Gtp3ecx/FadBXIQDM7j0Tn9fgiyBihDLWwrbwWAIuKY/iFxKmvr0d7e0eJIymlqTZ5HEM0dINQdcrUIhn45jkA9CkDOO2JFS+eZeyMYEyO0gboTJshMdhRdIjKUVz8+dnBMj0Y6BxByWE7AK1WeT1KGjpIZUIhiRUOh3VJbF5Xqn7pHJajOhbAlcXYAFgycElYkk26Gzj1SzqBuzMfDoa0VoXsdGlQy6AyO4CqTQaE5Tmo2uQAQIZMpAqQ6SNygkS8D40OUTo1pXSQ16JzlipQdp7M4+pu2B9v90jHJ4+X5GI7KcVIVpl1wZeMRCch+WIkirQmlhz5EXxQMotBjnpkkh3QKdnoK+I5JMmoFuKH2zSi5RvNzuExNLblm8/rylwymahIaUHVS7+cHHww4VC2GegIILNoiayyTJeDvHc5KJAZICQWZ19zIHHEEEu+mRylxTtE5ehK5i5JhyrVgXSs8ryEJEz8A5YDAek3YnJebm5ul2lqMuQSfy/yIyUkpW8kEkFDQ4OuxcAXh0SQapP3Lvfj9ejPki4Yttnn81mKuPH3RJH2Nhbf7ng1JlUjDXJKFAmZvsxUEpn5ICUNVRNVDO0xucop33R2AuOEckUuoJNAHExIo53/M4ogfWUALEXbOAKVtqWMA8Yb8gAsz4Gqj9tp7w0aNEiXm6QZccTECgn5JlFlyNEfJUy8fSX3i/eqS0gVKDufZOWQnA5GDhRkHJJSimpRSiF5PQ4m4mfg0A6k64Kef07hYkBbSmR5TdkO7iOlKtDp5JWZHdJ9c8SU45ZEkXaOdBfEG8Qyw4FqgPPt6NuRQWgp4SQZeH5KKgaYSVYpIfmXBKXKoRRjB1Iiyfgl2y8XZuKKYQUFBcjKytLxRLlyGW00qkBeP97rLm04zoOkt57/M835iFmvUAaF49UHYC3qKu0lOhlpDEsDWxKJ+8mwkDwP1RmXKCGJ2A5J/Pi8d8YWpY0n3Ra8P6AzUzYjIwODBw9Gfn6+vk5GRoYO/Uh3BEksC9NS0lLiyYA0nb38nf4xPg/5XBNBWhMrHvHxQABa1XWX9UBEIhFtYMtRG2ANdAOdncxzx6sHOmKlU5SSjVJJTkGToy6qOvqXZBiGBeLy8/Ph9/t1rJAe83A4jLy8PAtZ2H7ed/zIVGZAMIYZ7zbhtr6Ec4BvALHY8dJojXc/UAXEOzzZAXKoLdWkTCGhWpAjKSnRZAkk6UOTKdEkFUNDRx99tJ64wd+ZKQp0etqzs7N1O7Kzs+H1etHY2Kgdmqy+w/ujNIof9fF58F75EgSDQZ33xegDU36kxJLrCR0MaU8s+QZ25y2X36U6k8a6rAIjVSrfaJmOIlNS5PQtggTuzmsvJdDgwYN1KrEkqdvt1gU6MjIy4PP59HFU3ZzzyOsWFBR0KQsuU695fsBaVZrBdtpocvQnn6U8T6JIe2J1Z5/05G+RLgDpT5JxMRKMud5AZwkjeTyJJUdNPE+8HcZt9OL7/X49vYt5ViQEwzHyHAzPUBqSRPSK044KhUKWMuCAteKhdMFIm4mSiESS+WNsO1+ARPGNcJBKvxWALkSTM6S5jR0u7Rv+BlhdEdyfEoXSgz4tqjsZquF5eBwlks/nQyAQ6BKKYUqNy9VR7pvhIs6K5upfPC9Tj6VLgjHIeHLLZ8GBClUcXxpZz0EppUeZciEDWz3vb731FmbOnImSkhI4HA4899xzlt+VUrj99tsxdOhQZGVlYfr06di6datln9raWsyZM0c/5Hnz5um5cocCkoQf6VQEus98jA8HUV3IsEa85KPNIt0A3MZzyFglg968TlZWFnJycroN6LLTaahLCcEMCGk/8RjpuFRK6Yp/EiSEfA7SaSrbDHQOIijBZFgoUfSZWM3NzRg/fnyPtRfuvfdePPjgg3jsscewYcMG5OTkYMaMGZaiq3PmzMHmzZvx6quv4p///CfeeustzJ8/v69NsXjIJbqbCCHfRh4bTxDZYRz9xXvrmeMu1QMlmRzmywwH2jvs9PjcMXak9J5TklA9cjk5qW4pIUkAOSKV8Un5TOQKaLSxeJ8kKrfLdRoZVE8UDtUXP338wQ4Hnn32WVx00UUAOt6YkpIS/OQnP8FNN90EoCM3u6ioCCtXrsTll1+OTz/9FMcddxzef/99TJw4EQDw0ksv4bvf/S6++uorlJSUHPS6jY2N8Pv9GDFihMVu4JsFdLoZAFgevBxBSv8QJQpTiXk/8WnNfOhut1tPyGBHs7qwHMXRR+R0OpGbm6uLfVDSyTilXE+HlW2ysrKQm5trcWTy3kgsEp33yAmy0s6iapdqsa6uDk1NTZZUIOmr4n0wAyMSiaCyshINDQ163eie0K821vbt21FdXW0pXuv3+1FeXq6L165btw6BQECTCgCmT58Op9OJDRs2dHvenorbxquGeEkgR1JAZ2pvvKNPZm2yo+PFvrSdZHoJEW/DST8Y0FmVkMY3t0vvOPdnzA7oGvwmseRxlD7SIcy8KrksC++Dz5NB7HhVyOcVDof1p7m52aJ1DoZ+JRYL1BYVFVm2y+K11dXVGDJkiOV3t9uNgoKCHgvc9lTcVoZheB6Cb6/0wkvHoRz9SdUSTxDpGJV2G68vwy8EJ1Tww8C0VMssL8RRqbSrpHSULwklVDQa1flU0qCWLwhLL3EqGhcQZ458MBjUki9+BhIhXxyZVpQI0mJUuGjRIjQ0NOhPVVUVgK6jQZmDDqDbByEfpCSNnCsYnyYipYL8Kx+8tK9ISBKBM4ulPUcJw1QXHud0OrUXnKSXUoVk5L7xdiHQmV7DhL76+nocOHAAtbW1aGpqQnNzc5d2xj8jKXHjXTqJoF+JxQK18UuYyOK1xcXF2Ldvn+X39vZ21NbW9ljg1uPx6Ap7/BDSTyPVXryDsjtCSLdDvMdaPkgZEJbpNfJYmd/EDuYQXs5SlhKS98Z2UCrJa8TnQ8k2MRtBtoWkkj43SXigsz4Xj5WjQc4skpKrO2l2MPQrsUaMGIHi4mJUVFTobY2NjdiwYYMuXjtlyhTU19dj48aNep/XX38dsVgM5eXlfb4mOzw+nCPfOqI7svEv30o+XHYuj5OdQukk1SP9T7KGAiUW5+pJv1M8eZig53a79Xo8lFSs4UC/mZytw1QdzqKhGpR2HjM2SDLp1JVpO1y8gNeKH/Xa6nlvamrCtm3b9Pft27fjww8/REFBAYYNG4Yf/ehHuOuuuzBmzBiMGDECt912G0pKSvTIcdy4cTj33HNx9dVX47HHHkNbWxtuvPFGXH755QmNCCVIJD5EKY2k30eGNGgTcYTFfYHOEAcfPDtb2nH8X0oqOlrlIEGqWRlukZM4aFcBnTn6/MRiMWRnZ0MppW0i1gNlLI+k4kLqQCeZ+bIw7MR2yuRBuXIqiczSmDy3tEXjbcze0GdiffDBBzjjjDP0dxb2nzt3LlauXImbb74Zzc3NmD9/Purr6zF16lS89NJLlnVY/vKXv+DGG2/EWWedBaezo9Dtgw8+2NemWGwsmQkKdEqWeAM3/njAWkJSvtF8sPE5WZKcsgyjbAcJLtNxpNphR8b70Tj9SmaikrgkCtWrPK90eMqpZvJe5YsnJTvbk5eXh5qami5SLN7mTKhvDsePlSzQjzV69GgAnUaoy+XSsbf4oLK0E/jQJJkIkpOqjXabdE5Kv5AkEaUOnZAkPuf6UXJKVSgD4xy50TdFXxkdlZzxnJHRseJ8VlaW9ovRCOe6OLyn+HQhKdnli+hwdGRRsGwSl/GVKTexWAxffvllQn6stA5CS9LQ9pAzWYDe62PJ4DDtG/nAw+GwrhRDcP94VUpJx9RinqM79QxYKxHTjmIuO68hMx+kOpYOTZmNINUpIdspkxblOtFsJ6UWl66Tbpj4xMWDIa2JRRLJUZAU7dJYjzfmqZZkZwOwvM1USyQADWYpqWSAmvaNHBlKgvOatLekM5R1saTN5vF4tCTmgIKdzSJqTqdTr2pRV1enXQk8B++RWRNUo5KkJCjvhXagvG85nzERpDWxZIxL6n851I5PIYkP4sYb8NwubRZ2CGCtCcHOYcfQvmEHcegfiUSQl5enfVpyNErPtpSqjPUBVmNckpWjT26nE5RtlAVsOY2L8yBJboazKMllyrasCcaXKGme92RAeqblX2k6UirJcI0cDXV3Thk6obdahnN4XqnqGKCWGaZUKSQMXwS6NqT9IwPI8TOLZAipO3uNnnS2g+dj/VB+5MvFdsl0a2mbyWcr7zsRpLXEklJKGtjdxfoAWDpdqsD49Bn+pcQjORyOzlwmOZOHviNeg5JKOhxZHZkdKv1m0qXB6VqyLdKmAjoLqvFvKBTSNd6ZAEgbStpVshqhDN7Tkx/vlpHqktsSRVoTC4BFetDWkn6meEOe+8oRoVSJgDXTEoB2TDJEwo5n3XY5MJBqgxWOW1tbtTSRbgnWXpDSQroxZJtIGI54qQKDwSAOHDigCcHqftJBKtOVpd1JqRQfnpJhLbY5Puh+0H5JeM8UBR+AfAu787fEh264X3ymqAzPAJ2qId5rTTUpnZ4ypCPVCG0ppuTQPUCjW+5LSSHnOrJt0lMPdDiruZwuJZSsd0oJKP1c8aYCJZ/8nWSTI2Cq1USR1jZWfJK/JJQkT7xxH58qwofMTpQpKfI8csYLJQLVo7SbZP4VXRBM5pOecNo8gLWaiwz7ANAzZ3iN9vaOuly1tbV6cOHz+ZCXl6cluPS7SfUn7TlJUl6P9p1008Q7YRNB2kusnkaGEiSFfOvjHYjxEq07sS9Het2l6rAzeC75hnN9HKmGuYQK6y9IXxzQOd2eI1C2LxQKaWM9EAhYfpPqnv9LtSptSCm1gK7T6KmS+Xz64sdKa4kFWIPQ0qYgMdjx8Ql/QCdR4lVkdw8wvtNlCIbnz8jI0FmXVJUMZTHep5SyBI1pN1F6ycwJdibvgVVm6urqEIvFdNBbSlqgszw33Q5SesXn/0vnKUeH8STisUeMxOIDiDfO4zM7pStCxr7iY43yfExzlpkNVGsE94nPC5dBcWm3sbQ2R4eUoFICykhCfP45U5/liI8vFl0dvC9KJLlGEH1aUv1TKnE0KudJxqu/vszSSWtiJSKaZdBVDpfpeiCkSuTv0tEqbRPuK+N8QKfqig+/cqTX3NyMAwcOWDIXpH3H89JOk9Po2Sa5OICUUGwfbSraS/Ezc0jWeEcxv0sy8X5lvlaiSEtiyc6QIro7A13aEtyXvinAmqLMhyhtJZ5L+o0Y4pCGPa8TX6mGUo7nZa4TV9SiV5wEii/hLR2m8W6B+IEGJbAMKXGgwH34kvQ06qW/K34/KeUSyVtIS2JxebTPP/88yS05MhEMBuH3+3vdJy2JVVBQAADYtWvXQW/wSEZjYyNKS0tRVVV10DSXRKCUQjAYTCghMy2JRRXF5dgMekf8PIHDQaIvctq7GwxSE4ZYBrYgLYnl8XiwePFiy2oLBl2RzOeUljnvBqmPtJRYBqkPQywDW2CIZWALDLEMbEFaEmv58uUoKyuD1+tFeXk53nvvvWQ3acCQiqU6u0PaEevpp5/GwoULsXjxYmzatAnjx4/HjBkzulSw+aYilUp19gqVZpg0aZK64YYb9PdoNKpKSkrU0qVLk9iq5ACAevbZZ/X3WCymiouL1X333ae31dfXK4/Ho5588kmllFKffPKJAqDef/99vc+aNWuUw+FQu3fv7re2pZXEikQi2Lhxo6UUpdPpxPTp03UpyiMZdpXqPBSkFbEOHDiAaDTaaynKIxl2leo8FKQVsQzSB2lFrMGDB8PlcvVaivJIhl2lOg8FaUWszMxMTJgwwVKKMhaLoaKiQpeiPJKRjFKdPaLfhgEDhKeeekp5PB61cuVK9cknn6j58+erQCCgqqurk920AUEwGFSVlZWqsrJSAVD333+/qqysVDt37lRKKXX33XerQCCgnn/+efXxxx+r733ve2rEiBGqpaVFn+Pcc89V3/rWt9SGDRvU22+/rcaMGaNmz57dr+1MO2IppdRDDz2khg0bpjIzM9WkSZPU+vXrk92kAcPatWsVgC6fuXPnKqU6XA633XabKioqUh6PR5111llqy5YtlnPU1NSo2bNnq9zcXOXz+dRVV12lgsFgv7bTpM0Y2IK0srEM0geGWAa2wBDLwBYYYhnYAkMsA1tgiGVgCwyxDGyBIZaBLTDEMrAFhlgGtsAQy8AWGGIZ2AJDLANbYIhlYAsMsQxsgSGWgS0wxDKwBYZYfcDKlSvhcDiwY8eOPh03bdo0nHDCCf3alrKyMlx55ZX9es7+hCGWgS0wxDKwBYZYBrbAEOsw8Pzzz+P8889HSUkJPB4PRo0ahZ///Oc9LnG7ceNGnHrqqcjKysKIESPw2GOPddknHA5j8eLFGD16NDweD0pLS3HzzTcjHA7bfTv9irRcmSJVsHLlSuTm5mLhwoXIzc3F66+/jttvvx2NjY247777LPvW1dXhu9/9Li699FLMnj0bq1atwnXXXYfMzEz84Ac/ANAxq/vCCy/E22+/jfnz52PcuHH497//jWXLluHzzz/vUmQtpdGvsxS/4XjiiScUALV9+3allFKhUKjLPtdcc43Kzs5Wra2tetvpp5+uAKhf//rXels4HFYnn3yyGjJkiIpEIkoppf70pz8pp9Op/vWvf1nO+dhjjykA6p133tHbhg8frieppiKMKjwMZGVl6f+5mvxpp52GUCiEzz77zLKv2+3GNddco79nZmbimmuuwb59+3QdhdWrV2PcuHE49thjceDAAf0588wzAQBr164dgLvqHxhVeBjYvHkzfvazn+H1119HY2Oj5beGhgbL95KSEuTk5Fi2jR07FgCwY8cOTJ48GVu3bsWnn36KwsLCbq+XTuUwDbEOEfX19Tj99NPh8/lw5513YtSoUfB6vdi0aRNuueWWPq2fTMRiMZx44om4//77u/29tLT0cJs9YDDEOkS88cYbqKmpwd///nd85zvf0du3b9/e7f579uzRK9kTXMizrKwMADBq1Ch89NFHOOuss/q0YnwqwthYhwi5mj0RiUTwyCOPdLt/e3s7VqxYYdl3xYoVKCwsxIQJEwAAl156KXbv3o3f/va3XY5vaWlBc3Nzf96CrTAS6xBx6qmnIj8/H3PnzsUPf/hDOBwO/OlPf+pxveSSkhLcc8892LFjB8aOHYunn34aH374IR5//HFkZGQAAK644gqsWrUK1157LdauXYtvf/vbiEaj+Oyzz7Bq1Sq8/PLLlqK0KY1kD0vTCfHuhnfeeUdNnjxZZWVlqZKSEnXzzTerl19+WQFQa9eu1cedfvrp6vjjj1cffPCBmjJlivJ6vWr48OHq4Ycf7nKNSCSi7rnnHnX88ccrj8ej8vPz1YQJE9SSJUtUQ0OD3i/V3Q2mPpaBLTA2loEtMMQysAWGWAa2IKnEOpJX8fqmI2nEOtJX8fqmI2mjwvLycpxyyil4+OGHAXSEM0pLS7FgwQLceuutyWiSQT8iKQ5SruK1aNEiva0vq3jFYjHs2bMHeXl5aR/6SCcopRAMBlFSUgKns3dllxRi9baKV3y6CdCRVSkzKHfv3o3jjjvO9nYadI+qqiocffTRve6TFqPCpUuXwu/3648hVXKRl5d30H2SQqy+ruK1aNEiNDQ06E9VVdVANdWgGyRifiSFWH1dxcvj8cDn81k+BimOZAUpD2cVr4aGhm4XKjKfgfnIYHhPSGp2w6Gu4mWIlfrESsvshsbGRvj9/mQ344hFQ0PDQc2RtBgVGqQfDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkMsA1tgiGVgCwyxDGyBIZaBLTDEMrAFhlgGtsAQy8AWGGIZ2AJDLANbYIhlYAsMsQxsgSGWgS0wxDKwBYZYBrbAEMvAFhhiGdgCU+c9RVBWVga3u/vuqK2tRW1t7QC36DBxeHOZk4Nv2kzoadOmqZqamh7v9x//+IeaP3++cjgcSW8rEpwJbYiV5M/UqVPVjh07DnrPbW1t6h//+Ie64IILkt7mpBBr8eLFXRpyzDHH6N9bWlrU9ddfrwoKClROTo665JJLEioEIpGqxCotLVUTJ05Uo0aNSmj/yZMnq2Aw2Kd7b2lpUeedd96RSazjjz9eff311/qzf/9+/fu1116rSktLVUVFhfrggw/U5MmT1amnntqna6QisUaOHKlXQP3444/VySef3Ov+Z5xxhvr6668P6Rk3NDSomTNnHnnEGj9+fLe/1dfXq4yMDLV69Wq97dNPP1UA1Lp16xK+RioS68ILL7S08fPPP1dDhgzpcf977rnnkJ4vsXfvXnX++eenLLFscTds3boVJSUlGDlyJObMmYNdu3YB6FjFva2tDdOnT9f7HnvssRg2bFivRW3D4TAaGxstn1RCdnY2nnzyScu2MWPGwOPx2HbNIUOG4LnnnsOQIUNsu8bhoN+JVV5ejpUrV+Kll17Co48+iu3bt+O0005DMBhEdXU1MjMzEQgELMcUFRWhurq6x3PG1yBNtZVG58+fD6/XO+DXdbvdKVs1ut/9WOedd57+/6STTkJ5eTmGDx+OVatWWRbn7gsWLVqEhQsX6u+NjY0pQS632438/HxcfvnlBy1PbRcCgUCXWq6pANufRiAQwNixY7Ft2zYUFxcjEomgvr7esk9PRW2JVKxBOnXqVNxyyy3Yu3cvysvLu/z+3nvvIRQK2d6OZ5991vZrHAps97w3NTXhiy++wBVXXIEJEyYgIyMDFRUVmDVrFgBgy5Yt2LVrV7dFbVMN8+bNw6mnngoAOPvss3uVmo899hhqampsb9MRowpvuukmzJw5E8OHD8eePXuwePFiuFwuzJ49G36/H/PmzcPChQtRUFAAn8+HBQsWYMqUKZg8eXJ/N+WQUFZWptVaeXk57rnnHv3boEGDkJ2dfdBz/OMf/8CqVatsa2M6oN+J9dVXX2H27NmoqalBYWEhpk6divXr16OwsBAAsGzZMjidTsyaNQvhcBgzZszocYHugcbMmTOxevXqwx7NhUKhtFoY3BYcljMlSbDLj/XCCy8opZSKxWIqGo3qT01Njbr44ovVxRdfrA4cOKC394Sqqip11lln9Xqtw/VjEZs3b05JP5bJbhD4/ve/jxdffBFPP/00/vjHP+rtSik0NDQAAF5//XU4HA64XC4899xzGDVqFIYOHWo5z9FHH41hw4YNSJsvvPDCAblOn9Evr80AI5U872eccYaqqqrq0sZnnnlG5eXl2S6xioqKUlJimUS/w8TatWvxH//xH11cC7NmzcLKlSttvXYkEoFK0TL9hlj9gPXr1+Oyyy7rsv3kk0+29brz589P2RVpDbH6CZFIpMu2FStW9Lj/iy++iK+//vqQr7dx40Zs2LDhkI+3G4ZY/YS2tja0trYC6CDZokWLsGzZsh73f/PNN7Ft27ZDulZDQwPWrVvX7aKhKYN+sSAHGKlkvMvPlVdeqfbt26d+/OMfJ7T/oEGDenVbdIfKykpVUlKS1DRlk5qchE9fUoc9Ho+6+uqr1ZtvvpnQfX/00UfquOOOS/o9mtW/0gSDBg3CCy+8gKFDh2Lw4MHIzs5GMBhEXV0dAOB3v/sd/vSnPyEUCqWEsZ7I6l+GWCmGa665BhMnTsQrr7yC1atXJ7s53cIQy8AWmPUKDZIGQywDW2CIZWALDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkMsA1tgiGVgCwyxDGyBIZaBLTDEMrAFhlgGtqDPxHrrrbcwc+ZMlJSUwOFw4LnnnrP8rpTC7bffjqFDhyIrKwvTp0/H1q1bLfvU1tZizpw58Pl8CAQCmDdvHpqamg7rRgxSC30mVnNzM8aPH4/ly5d3+/u9996LBx98EI899hg2bNiAnJwczJgxQ89gAYA5c+Zg8+bNePXVV/HPf/4Tb731FubPn3/od2GQejjkGQ0dmafq2Wef1d9jsZgqLi5W9913n95WX1+vPB6PevLJJ5VSSn3yyScKgHr//ff1PmvWrFEOh0Pt3r07oeum8mSKI+Ez4FPst2/fjurqakvxWr/fj/Lycl28dt26dQgEApg4caLeZ/r06XA6nT1OwEz14rYGXdGvxGKB2qKiIst2Wby2urq6S6Vft9uNgoKCHgvcpnpxW4OuSItR4aJFi9DQ0KA/VVVVyW6SwUHQr8Rigdr4Kr6yeG1xcXGXuXHt7e2ora3tscBtKha3Negd/UqsESNGoLi4GBUVFXpbY2MjNmzYoIvXTpkyBfX19di4caPe5/XXX0csFuu2+rBBmiLBAaBGMBhUlZWVqrKyUgFQ999/v6qsrFQ7d+5USil19913q0AgoJ5//nn18ccfq+9973tqxIgRqqWlRZ/j3HPPVd/61rfUhg0b1Ntvv63GjBmjZs+enXAbzKgw9UeFfSbW2rVru73Y3LlzlVIdLofbbrtNFRUVKY/Ho8466yy1ZcsWyzlqamrU7NmzVW5urvL5fOqqq67q0ypYhlipTywzE9qgzzAzoQ2SBkMsA1tgiGVgCwyxDGyBIZaBLTDEMrAFhlgGtsAQy8AWGGIZ2AJDLANbYIhlYAsMsQxsQVoSKw3j5t8oJPL805JYA7E6vEHPCAaDB90nLZfuLSgoAADs2rXLpM/0gsbGRpSWlqKqqqpf0rmVUggGgygpKTnovmlJLKezQ9D6/X6T/54A+nOeQKIvclqqQoPUhyGWgS1IS2J5PB4sXrwYHo8n2U1JaSTzOaVlzrtB6iMtJZZB6sMQy8AWGGIZ2AJDLANbkJbEWr58OcrKyuD1elFeXo733nsv2U0aMKRLqc60I9bTTz+NhQsXYvHixdi0aRPGjx+PGTNmpMTq7gOBtCnVmXDBhBTBpEmT1A033KC/R6NRVVJSopYuXZrEViUHQHJKdSaCtJJYkUgEGzdutJSidDqdmD59ui5FeSTDrlKdh4K0ItaBAwcQjUZ7LUV5JMOuUp2HgrQilkH6IK2INXjwYLhcrl5LUR7JsKtU56EgrYiVmZmJCRMmWEpRxmIxVFRU6FKURzJSqlRnvw0DBghPPfWU8ng8auXKleqTTz5R8+fPV4FAQFVXVye7aQOCVCjVmQjSjlhKKfXQQw+pYcOGqczMTDVp0iS1fv36ZDdpwJAKpToTgUmbMbAFaWVjGaQPDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkMsA1tgiGVgCwyxDGyBIZaBLTDEMrAFhlgGtsAQy8AWGGIZ2AJDLANbYIhlYAsMsQxsgSGWzXjkkUewcuXKAbnWJ598gjvuuAM7duwYkOv1BkMsmzHQxFqyZIkhlsE3GP065+cbgk2bNqlzzz1X5eXlqZycHHXmmWeqdevW6d8XL16sunt0TzzxhAKgtm/frpRSavjw4V2maZ1++umWfd988001f/58VVBQoPLy8tQVV1yhamtrLecFoBYvXtzlesOHD9fTvni++M/atWv745H0GWm5MoWd2Lx5M0477TT4fD7cfPPNyMjIwIoVKzBt2jS8+eabfZot/MADD2DBggXIzc3F//t//w9A14IdN954IwKBAO644w5s2bIFjz76KHbu3Ik33ngDDocj4Wt95zvfwQ9/+EM8+OCD+J//+R+MGzcOAPTfAUdS6JzCuOiii1RmZqb64osv9LY9e/aovLw89Z3vfEcplbjEUkqp448/Xkup7vadMGGCikQievu9996rAKjnn39eb0MCEksppVavXp1UKSVhbCyBaDSKV155BRdddBFGjhyptw8dOhTf//738fbbb6OxsbFfrzl//nxkZGTo79dddx3cbjdefPHFfr3OQMMQS2D//v0IhUI45phjuvw2btw4xGIxVFVV9es1x4wZY/mem5uLoUOHpsTI7nBgiHUI6Mn2iUajA9qOgb5eX2CIJVBYWIjs7Gxs2bKly2+fffYZnE4nSktLkZ+fDwCor6+37LNz584uxx3MAI+vaNzU1ISvv/4aZWVlelt+fn6Xa0UiEXz99dd9utZAwhBLwOVy4ZxzzsHzzz9vUUV79+7FX//6V0ydOhU+nw+jRo0C0FEam2hubsYf/vCHLufMycnpQgqJxx9/HG1tbfr7o48+ivb2dpx33nl626hRoyzX4nHxEisnJwdAV8InA8bdEIe77roLr776KqZOnYrrr78ebrcbK1asQDgcxr333gsAOOecczBs2DDMmzcPP/3pT+FyufC///u/KCwsxK5duyznmzBhAh599FHcddddGD16NIYMGYIzzzxT/x6JRHDWWWfh0ksvxZYtW/DII49g6tSpuPDCC/U+//Vf/4Vrr70Ws2bNwtlnn42PPvoIL7/8MgYPHmy51sknnwyXy4V77rkHDQ0N8Hg8OPPMM7vUHB0QJHtYmorYtGmTmjFjhsrNzVXZ2dnqjDPOUO+++65ln40bN6ry8nKVmZmphg0bpu6///5u3Q3V1dXq/PPPV3l5eb06SPPz81Vubq6aM2eOqqmpsVwrGo2qW265RQ0ePFhlZ2erGTNmqG3btnVxNyil1G9/+1s1cuRI5XK5kup6MPWxkoSVK1fiqquuwvvvv28pjf1NgbGxDGyBIZaBLTDEMrAFSSXWkbyK15VXXgml1DfSvgKSSKwjfRWvbzqSNiosLy/HKaecgocffhhAx0IApaWlWLBgAW699dZkNMmgH5EUBylX8Vq0aJHe1pdVvGKxGPbs2YO8vLyUCmN806GUQjAYRElJCZzO3pVdUojV2ypen332WZf9w+EwwuGw/r57924cd9xxtrfToHtUVVXh6KOP7nWftBgVLl26FH6/X38MqZKLvLy8g+6TFGL1dRWvRYsWoaGhQX/6OyfKoG9IxPxICrH6uoqXx+OBz+ezfAxSHEmJUKrDW8WroaGh2xkp5jMwn4aGhoP2UVKzGw51FS9DrNQnVlpmNzQ2NsLv9ye7GUcsGhoaDmqOpMWo0CD9YIhlYAsMsQxsgSGWgS0wxDKwBYZYBhrZ2dm49NJL++VchlgGGpFIpN+SLQ2xDDTa29v7rWaEIZaBLTDEMugVS5YsQW5ubp+PMyEdgz7DhHQMkgZDLINuMXz48MM63hDLoFtccMEFh3W8sbEM+gxjYxkkDYZYBrbAEMtAIysrq9/OZYhlAAA49thjsXTp0n47nyGWAYCOuZ533nlnv53PFLc9wuFwOOD3+/H222/363mNxDrCkZWVhddeew3HHntsv57XEOsIRzgcxmuvvWapK98fMA7SFEFZWRncbjeWLFmCYcOG4aqrrgIAtLa24quvvrL12h6PB263G83NzQntn4iDNC3rvH/TZkJPmzatS213Ytu2bWrixIlJb6P8pPwU+0PFN4VYubm56vnnn1c7duzo9X4//fRTNXbs2KS3N6nE4iKR8nPMMcfo31taWtT111+vCgoKVE5OjrrkkksSKgQikYrE8vl8auLEierkk09O+JhTTz1VxWKxhO5579696phjjkmJtieNWMcff7z6+uuv9Wf//v3692uvvVaVlpaqiooK9cEHH6jJkyerU089tU/XSEVinXPOOUoppZqbm9XFF1+c0DEffvhhwvfc2tqq/vM//zMl2p40Yo0fP77b3+rr61VGRoZavXq13vbpp58qAJbFvA+GVCOW1+tVlZWVun379+9X559/fq/HZGRkqI8//rjHe3zyySfVv/71L9Xc3KwikYhqb29Xv/rVr5Tb7e7Xtns8HkuVnz/84Q+pS6zs7Gw1dOhQNWLECPX9739f7dy5UymlVEVFhQKg6urqLMdwkaOe0NraqhoaGvSnqqoq6WSSn+zsbMu6zkop9etf/1p5vd4ej/nlL3/ZoxoMBoPq6quvVpMmTVJer1fNmzdPVVRUqGg0qn7zm9+orKysXtvjcrnU2Wefrc455xx1zjnnqBNOOKHX/cvLy9XLL7+sVq1alRBxk0KsF198Ua1atUp99NFH6qWXXlJTpkxRw4YNU42Njeovf/mLyszM7HLMKaecom6++eYez9md3ZZKnx/96EcqGo12aXdpaWmPx9xzzz093u/+/fvVzJkz1UUXXaQ72u/3q/b2dvW3v/1NBQKBXtvzk5/8xELaTz75pF9HlikxKqyrq1M+n0/97ne/O2RiparEcrvdqrCwsMeCcb0R67bbblPt7e093vOvfvUry/5Op1Mdf/zxqqCgoNc2XXfddSoUCnU53w033KCcTueAEct2z3sgEMDYsWOxbds2FBcXIxKJdFkBtKeitkSq1SAdPXo0LrvsMtxyyy3Yu3cvysvLu+zz3nvvIRQK9XiOn//85/i///u/bn979tln8eGHH1q2xWIxbN68GbW1tb22rbi4uNv0l4ceeqhL+XM7YXsQuqmpCV988QWuuOIKTJgwARkZGaioqMCsWbMAAFu2bMGuXbu6LWqbSpg8eTKuvvpqAMBxxx2HyZMn97jve++9hzlz5qCmpuaQrvX555/jz3/+8yEd2xMGeqGFfifWTTfdhJkzZ2L48OHYs2cPFi9eDJfLhdmzZ8Pv92PevHlYuHAhCgoK4PP5sGDBAkyZMqXXjhpIlJWVWVZdWLBgAWbNmoWcnBwUFBSgra0NVVVVCIVCyM7O7vYc27dvx7Zt2w7p+kopvPzyy4d0LADU1tYiHA7D4/FYtv/4xz/uUv7cTvQ7sb766ivMnj0bNTU1KCwsxNSpU7F+/XoUFhYCAJYtWwan04lZs2YhHA5jxowZeOSRR/q7GYeEmTNnYvXq1V06RaKurg733Xcfrr/+epx44on9du2mpiZs3LgRp59+On7zm9/gpJNOOqTz/OY3v0FmZiZGjhxp2b5u3TrEYrH+aGpiOKgVloKwy4/1wgsvKKU61mDm57e//a26+OKL9eeGG27Qv/WEp556KqHrSQdpKBRSmzZt0vc3d+7cLvtPnDhRzZgxI+mDlpQYFdoBu4jl9/vV448/rgKBgP54PJ4uo7NAIKAGDRqk/vWvf6k9e/Z0ad/LL7+sBg0a1CdixeOPf/yjXqCcn4yMDDV69OgubTLE6iekkuf9jDPOUFVVVV3a+Mwzz3QhhvxMmzatW1JKrFixoovDcvny5Wr06NGGWHYglYgFQE2ePFk1Nzd3aeff/va3Ho8ZOXKk2rZtW6/3GQ6HVVZWlnK73WrQoEHqlVdeUTt37kwLYpkM0n7A+vXrcdlll3XZ3pux/OWXX6KpqanX8zocDkybNg2///3vsXv3boTDYTzwwAO9+sdSBWYyRT8hEol02bZw4cLDOmdGRgZefPFF/T0Wi2HZsmWHdc6BgiGWTbjjjjtQXV3d6z6NjY0JnaulpQVXXnllei3GfsiGThKRajYWADV+/Hj1zDPPqH379qlnnnlGnXHGGQc9ZtCgQb26LZTqCEjPmjUr6fcnP8Z4T8LnggsuSHhfj8ejHn300R7vs62tLeVIZYiVJp/s7GxVVlam1q9fb5lQ0djYqGbNmqUcDkfS23goxDLTv1IIU6dOxdy5cwEAr7zyClavXp3kFnWPRKZ/GWIZ9Bmm8JpB0mCIZWALDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkMsA1tgiGVgCwyx0hA/+9nPkt2Eg8IQK81QXFyM7373u8luxkFhiJVmuO2227Bhw4ZkN+OgMKnJaYbly5f320rzdsJIrDRDKBQa2Knyh4g+E+utt97CzJkzUVJSAofDgeeee87yu1IKt99+O4YOHYqsrCxMnz4dW7dutexTW1uLOXPmwOfzIRAIYN68eQedCmXQgXHjxvXrKl12oc/Eam5uxvjx47F8+fJuf7/33nvx4IMP4rHHHsOGDRuQk5ODGTNmoLW1Ve8zZ84cbN68Ga+++ir++c9/4q233sL8+fMP/S6OIKxZswZ1dXXJbsbBcTi55wDUs88+q7/HYjFVXFys7rvvPr2tvr5eeTwe9eSTTyqlOsoWAlDvv/++3mfNmjXK4XCo3bt3J3Tdb1LOezp+Bnwm9Pbt21FdXY3p06frbX6/H+Xl5Vi3bh2AjnI6gUAAEydO1PtMnz4dTqezx9FOOBxGY2Oj5WOQ2uhXYnGCZnxJwqKiIv1bdXU1hgwZYvnd7XajoKCgxwmeS5cuhd/v15/S0tL+bLaBDUiLUeGiRYvQ0NCgP1VVVcluksFB0K/EYoHa+JKEsnhtcXEx9u3bZ/m9vb0dtbW1PRa4TbXitgYHR78Sa8SIESguLkZFRYXe1tjYiA0bNujitVOmTEF9fT02btyo93n99dcRi8W6rT5skKZIcACoEQwGVWVlpaqsrFQA1P33368qKyv16hN33323CgQC6vnnn1cff/yx+t73vqdGjBihWlpa9DnOPfdc9a1vfUtt2LBBvf3222rMmDFq9uzZCbfBjApTf1TYZ2KtXbu224vNnTtXKdXhcrjttttUUVGR8ng86qyzzlJbtmyxnKOmpkbNnj1b5ebmKp/Pp6666ioVDAYTboMhVuoTy8yENugzzExog6TBEMvAFhhiGdgCQywDW2CIZWALDLEMbIEhloEtMMQysAWGWAa2wBDLwBYYYhnYAkMsA1uQlsRKw7j5NwqJPP+0JNahrg5v0D8IBoMH3Sctp9gXFBQAAHbt2mXSZ3pBY2MjSktLUVVV1S/p3EopBINBlJSUHHTftCSW09khaP1+v8l/TwD9OU8g0Rc5LVWhQerDEMvAFqQlsTweDxYvXgyPx5PspqQ0kvmc0jLn3SD1kZYSyyD1YYhlYAsMsQxsgSGWgS1IS2ItX74cZWVl8Hq9KC8vx3vvvZfsJg0Y0qVUZ9oR6+mnn8bChQuxePFibNq0CePHj8eMGTO6VLD5piJtSnUmXDAhRTBp0iR1ww036O/RaFSVlJSopUuXJrFVyQGQnFKdiSCtJFYkEsHGjRstpSidTiemT5+uS1EeybCrVOehIK2IdeDAAUSj0V5LUR7JsKtU56EgrYhlkD5IK2INHjwYLper11KURzLsKtV5KEgrYmVmZmLChAmWUpSxWAwVFRW6FOWRjJQq1dlvw4ABwlNPPaU8Ho9auXKl+uSTT9T8+fNVIBBQ1dXVyW7agCAVSnUmgrQjllJKPfTQQ2rYsGEqMzNTTZo0Sa1fvz7ZTRowpEKpzkRg0mYMbEFa2VgG6QNDLANbYIhlYAsMsQxsgSGWgS0wxDKwBYZYBrbAEMvAFhhiGdgCQywDW2CIZWALDLEMbMH/B0prJ0RdeoG9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best IoU: 0.4327, Best Dice: 0.5713, Best Acc: 0.9524\n"
     ]
    }
   ],
   "source": [
    "# -----------------------------\n",
    "# Threshold method application\n",
    "# -----------------------------\n",
    "\n",
    "# empty the cache for model and training process\n",
    "torch.cuda.empty_cache()\n",
    "\n",
    "# Do the segmentation to exteract the spine part\n",
    "# Threshold segmentation hyper parameters\n",
    "minThreshold = 250\n",
    "maxThreshold = 1500\n",
    "\n",
    "# the variables for evaluation\n",
    "total_correct_pixels = 0\n",
    "total_pixels = 0\n",
    "total_iou = 0.0\n",
    "total_dice = 0.0\n",
    "\n",
    "# set the hyper parameters and the paths 设置超参数和路径\n",
    "data_path = \"./dataset/Spineweb_dataset15/train_9/\"    # the path of the dataset\n",
    "file_list = os.listdir(data_path)  # list all files and folders in the folder(data path)\n",
    "file_count = len(file_list)  # get the numer of files\n",
    "print(\"Length of the dataset:\",file_count)\n",
    "resize_scale = (128,128)    # set the resize scale for images\n",
    "print(\"----------------  Loading the dataset and evaluate  ----------------\")\n",
    "data_count = 0\n",
    "poster_threshold = 200\n",
    "for index in range(file_count):     # iterate the whole dataset\n",
    "    metadata = np.load(data_path+str(index)+\".npy\")     # load the data\n",
    "    image = np.array(Image.fromarray(metadata[0]).resize(resize_scale))     # Using the PIL to resize the transfer into npy again\n",
    "    label = Image.fromarray((metadata[1] * 255).astype('uint8')).resize(resize_scale)     # Using the PIL to resize the transfer into npy again\n",
    "    label = np.array(label)\n",
    "    label[label < poster_threshold] = 0\n",
    "    label[label >= poster_threshold] = 1\n",
    "    \n",
    "    # print(type(image),image.shape,image)\n",
    "    # print(type(label),label.shape,label)\n",
    "    output = CTimageLIB.thresholdSegmentation(image,minThreshold,maxThreshold)\n",
    "    # evaluations\n",
    "    current_correct_pixels,current_pixels = calculate_pixel_accuracy(output,label)\n",
    "    current_iou = calculate_iou(output,label)\n",
    "    current_dice = calculate_dice(output,label)\n",
    "    # add the results to total one\n",
    "    total_correct_pixels += current_correct_pixels\n",
    "    total_pixels += current_pixels\n",
    "    total_iou += current_iou\n",
    "    total_dice += current_dice\n",
    "    data_count += 1\n",
    "    # visualize the results\n",
    "    if index == 100:\n",
    "        visualize_results(image,label,output)\n",
    "\n",
    "final_acc = total_correct_pixels/total_pixels\n",
    "final_iou = total_iou/data_count\n",
    "final_dice = total_dice/data_count\n",
    "\n",
    "print(f\"Best IoU: {final_iou:.4f}, Best Dice: {final_dice:.4f}, Best Acc: {final_acc:.4f}\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dataEngineering",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}