tumor_no_tumor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8f5e879fc9b545819a0d4ff00d5ad213",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "0it [00:00, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GLCM is symmetrical, therefore Sum Average = 2 * Joint Average, only 1 needs to be calculated\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "129\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.ndimage import label, find_objects\n",
    "import matplotlib.pyplot as plt\n",
    "from pathlib import Path\n",
    "from tqdm.auto import tqdm\n",
    "import SimpleITK as sitk\n",
    "import pandas as pd\n",
    "import radiomics\n",
    "from radiomics import featureextractor\n",
    "import os\n",
    "\n",
    "def count_blob_pixels(array):\n",
    "    # Label the blobs in the array\n",
    "    labeled_array, num_features = label(array)\n",
    "\n",
    "    # Calculate the number of pixels in each blob\n",
    "    blob_sizes = [(labeled_array == i).sum() for i in range(1, num_features + 1)]\n",
    "\n",
    "    return blob_sizes\n",
    "\n",
    "\n",
    "def remove_small_blobs(array, min_size=100):\n",
    "    # Label the blobs in the array\n",
    "    labeled_array, num_features = label(array)\n",
    "\n",
    "    # Find the sizes of each blob\n",
    "    blob_sizes = [(labeled_array == i + 1).sum() for i in range(num_features)]\n",
    "\n",
    "    # Create a mask to keep blobs that meet the size requirement or if it's the only blob\n",
    "    if num_features == 1:\n",
    "        mask = labeled_array != 0  # Keep the single blob regardless of size\n",
    "    else:\n",
    "        mask = np.isin(\n",
    "            labeled_array,\n",
    "            [i + 1 for i, size in enumerate(blob_sizes) if size >= min_size],\n",
    "        )\n",
    "\n",
    "    # Apply the mask to the original array\n",
    "    filtered_array = array * mask\n",
    "\n",
    "    return filtered_array\n",
    "\n",
    "\n",
    "def extract_radiomics_features(pet_img, mask_array, output_dir, sub_name, label_type):\n",
    "    if not os.path.exists(output_dir):\n",
    "        os.makedirs(output_dir)\n",
    "    mask_array = remove_small_blobs(mask_array)\n",
    "    # Label the blobs in the array\n",
    "    labeled_array, num_features = label(mask_array)\n",
    "\n",
    "    # Initialize the PyRadiomics feature extractor\n",
    "    extractor = featureextractor.RadiomicsFeatureExtractor()\n",
    "\n",
    "    for i in range(1, num_features + 1):\n",
    "        # Create the filename\n",
    "\n",
    "        filename = f\"{sub_name}_blob_{i}_label_{label_type}.csv\"\n",
    "\n",
    "        if Path(os.path.join(output_dir, filename)).exists():\n",
    "            continue\n",
    "        blob_mask = (labeled_array == i).astype(np.uint8)\n",
    "        if blob_mask.sum()<=50:\n",
    "            continue\n",
    "        # Ensure the blob mask has the same shape as the original array\n",
    "        assert blob_mask.shape == mask_array.shape\n",
    "        mask_img = sitk.GetImageFromArray(blob_mask)\n",
    "        mask_img.CopyInformation(pet_img)\n",
    "        # Extract features\n",
    "        feature_vector = extractor.execute(pet_img, mask_img)\n",
    "        print(len(feature_vector))\n",
    "        # Convert the feature vector to a pandas DataFrame\n",
    "        feature_df = pd.DataFrame([feature_vector])\n",
    "\n",
    "        # Remove non-numeric columns\n",
    "        feature_df = feature_df.select_dtypes(include=[np.number])\n",
    "\n",
    "        # Save the DataFrame as a CSV file\n",
    "        feature_df.to_csv(os.path.join(output_dir, filename), index=False)\n",
    "\n",
    "\n",
    "def calculate_tp_fp_fn(gt, pred):\n",
    "    # Label the blobs\n",
    "    gt_labels, num_gt = label(gt)\n",
    "    pred_labels, num_pred = label(pred)\n",
    "\n",
    "    # Find bounding boxes of blobs\n",
    "    gt_slices = find_objects(gt_labels)\n",
    "    pred_slices = find_objects(pred_labels)\n",
    "\n",
    "    # Create empty masks for TP, FP, FN\n",
    "    tp_mask = np.zeros_like(gt)\n",
    "    fp_mask = np.zeros_like(gt)\n",
    "    fn_mask = np.zeros_like(gt)\n",
    "\n",
    "    # Calculate TP, FP, FN\n",
    "    tp, fp, fn = 0, 0, 0\n",
    "\n",
    "    # Helper function to determine if two blobs overlap\n",
    "    def overlap(slice1, slice2):\n",
    "        return not (\n",
    "            slice1[0].stop < slice2[0].start\n",
    "            or slice1[0].start > slice2[0].stop\n",
    "            or slice1[1].stop < slice2[1].start\n",
    "            or slice1[1].start > slice2[1].stop\n",
    "            or slice1[2].stop < slice2[2].start\n",
    "            or slice1[2].start > slice2[2].stop\n",
    "        )\n",
    "\n",
    "    matched_gt = set()\n",
    "    matched_pred = set()\n",
    "\n",
    "    for i, gt_slice in enumerate(gt_slices):\n",
    "        gt_blob = gt_labels == (i + 1)\n",
    "        match_found = False\n",
    "        for j, pred_slice in enumerate(pred_slices):\n",
    "            if overlap(gt_slice, pred_slice):\n",
    "                pred_blob = pred_labels == (j + 1)\n",
    "                tp_mask[gt_blob & pred_blob] = 1  # True positive overlap region\n",
    "                tp += 1\n",
    "                matched_gt.add(i)\n",
    "                matched_pred.add(j)\n",
    "                match_found = True\n",
    "                break\n",
    "        if not match_found:\n",
    "            fn_mask[gt_blob] = 1\n",
    "\n",
    "    for j, pred_slice in enumerate(pred_slices):\n",
    "        if j not in matched_pred:\n",
    "            pred_blob = pred_labels == (j + 1)\n",
    "            fp_mask[pred_blob] = 1\n",
    "\n",
    "    fn = num_gt - tp\n",
    "    fp = num_pred - tp\n",
    "\n",
    "    return tp_mask, fp_mask, fn_mask, tp, fp, fn\n",
    "\n",
    "\n",
    "src_dir = Path(\"/mnt/nfs/slow_ai_team/organ_segmentation/nnunet_liverv0.0/nnUNet_raw_database/nnUNet_raw/nnUNet_raw_data/Dataset019_AutoPET2024/\")\n",
    "pred_dir = \"imagesTr_pred_res\"\n",
    "gt_dir = \"imagesTr\"\n",
    "metrics=[]\n",
    "for pred in tqdm((src_dir/pred_dir).glob(\"*.nii.gz\")):\n",
    "    fname = (pred.name).split(\".nii\")[0]\n",
    "    gt = str(pred).replace(\"imagesTr_pred_res\", \"labelsTr\")\n",
    "    pred_array = sitk.GetArrayFromImage(sitk.ReadImage(pred))\n",
    "    pet_img = (sitk.ReadImage(\n",
    "        src_dir / \"imagesTr\" / str(fname + \"_0000.nii.gz\")\n",
    "    ))\n",
    "    gt_array = sitk.GetArrayFromImage(sitk.ReadImage(gt))\n",
    "    gt_array = gt_array==20\n",
    "    pred_array = pred_array==20\n",
    "    tp_mask, fp_mask, fn_mask, tp, fp, fn = calculate_tp_fp_fn(gt_array, pred_array)\n",
    "    fp_voxels = count_blob_pixels(fp_mask)\n",
    "    fn_voxels = count_blob_pixels(fn_mask)\n",
    "    tp_voxels = count_blob_pixels(tp_mask)\n",
    "    radiomics_path = \"tp_fp_fn_radiomics_ct\"\n",
    "    extract_radiomics_features(pet_img, fp_mask, radiomics_path, fname, \"fp\")\n",
    "    extract_radiomics_features(pet_img, tp_mask, radiomics_path, fname, \"tp\")\n",
    "    extract_radiomics_features(pet_img, fn_mask, radiomics_path, fname, \"fn\")\n",
    "    metrics.append({\"Subject\": pred.name, \"TP\":tp, \"FP\": fp, \"FN\": fn, \"TP_voxels\": tp_voxels, \"FP_voxels\": fp_voxels, \"FN_voxels\": fn_voxels})\n",
    "df = pd.DataFrame(metrics)\n",
    "df.head()\n",
    "df.to_csv(\"TP_FP_FN_res_model.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   diagnostics_Image-original_Mean  diagnostics_Image-original_Minimum  \\\n",
      "0                      -908.964342                        -1368.807739   \n",
      "1                      -912.161126                        -1457.807739   \n",
      "2                      -881.026459                        -1598.900391   \n",
      "3                      -957.705830                        -1313.069824   \n",
      "4                      -919.534544                        -1437.792358   \n",
      "\n",
      "   diagnostics_Image-original_Maximum  diagnostics_Mask-original_VoxelNum  \\\n",
      "0                         3265.816895                                 138   \n",
      "1                         3243.423096                                 570   \n",
      "2                         3737.096680                                 724   \n",
      "3                         3330.528320                                 189   \n",
      "4                         3217.746582                                 115   \n",
      "\n",
      "   diagnostics_Mask-original_VolumeNum  original_shape_Elongation  \\\n",
      "0                                    1                   0.912729   \n",
      "1                                    1                   0.411556   \n",
      "2                                    1                   0.800805   \n",
      "3                                    1                   0.793835   \n",
      "4                                    1                   0.800569   \n",
      "\n",
      "   original_shape_Flatness  original_shape_LeastAxisLength  \\\n",
      "0                 0.773735                       11.569314   \n",
      "1                 0.277550                       13.021757   \n",
      "2                 0.543631                       45.277825   \n",
      "3                 0.721216                       12.832862   \n",
      "4                 0.684007                       10.490139   \n",
      "\n",
      "   original_shape_MajorAxisLength  original_shape_MinorAxisLength  \\\n",
      "0                       14.952549                       13.647625   \n",
      "1                       46.916834                       19.308902   \n",
      "2                       83.287803                       66.697295   \n",
      "3                       17.793379                       14.125013   \n",
      "4                       15.336298                       12.277758   \n",
      "\n",
      "   original_shape_VoxelVolume  label  \n",
      "0                 1716.860837      1  \n",
      "1                 7091.381716      1  \n",
      "2                60048.363939      1  \n",
      "3                 2351.352885      1  \n",
      "4                 1430.717364      1  \n",
      "2000\n"
     ]
    }
   ],
   "source": [
    "import glob\n",
    "def load_and_label_radiomics_data(folder_path, num_files=1000):\n",
    "    # Get the list of tp and fp files\n",
    "    tp_files = glob.glob(os.path.join(folder_path, \"*_tp.csv\"))[:num_files]\n",
    "    fp_files = glob.glob(os.path.join(folder_path, \"*_fp.csv\"))[:num_files]\n",
    "\n",
    "    # Load tp files into a DataFrame and add the label column\n",
    "    tp_dataframes = []\n",
    "    for file in tp_files:\n",
    "        df = pd.read_csv(file)\n",
    "        df[\"label\"] = 1\n",
    "        tp_dataframes.append(df)\n",
    "\n",
    "    # Load fp files into a DataFrame and add the label column\n",
    "    fp_dataframes = []\n",
    "    for file in fp_files:\n",
    "        df = pd.read_csv(file)\n",
    "        df[\"label\"] = 0\n",
    "        fp_dataframes.append(df)\n",
    "\n",
    "    # Concatenate all DataFrames\n",
    "    all_data = pd.concat(tp_dataframes + fp_dataframes, ignore_index=True)\n",
    "\n",
    "    return all_data\n",
    "\n",
    "\n",
    "# Example usage\n",
    "# folder_path = \"tp_fp_fn_radiomics\"\n",
    "combined_data = load_and_label_radiomics_data(radiomics_path, num_files=1000)\n",
    "print(combined_data.head())\n",
    "print(len(combined_data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classification Report:\n",
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       0.71      0.05      0.09       201\n",
      "           1       0.51      0.98      0.67       199\n",
      "\n",
      "    accuracy                           0.51       400\n",
      "   macro avg       0.61      0.51      0.38       400\n",
      "weighted avg       0.61      0.51      0.38       400\n",
      "\n",
      "Accuracy: 0.51\n",
      "Confusion Matrix:\n",
      "[[ 10 191]\n",
      " [  4 195]]\n",
      "[0.10261634 0.09822058 0.1035442  0.12292938 0.12439349 0.12003546\n",
      " 0.11487313 0.10595383 0.10743359]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import ExtraTreesClassifier, RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import classification_report, accuracy_score, confusion_matrix\n",
    "from sklearn.feature_selection import VarianceThreshold\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "def classify(cdf):\n",
    "    y = cdf[\"label\"]\n",
    "    X = cdf.drop(columns=[\"label\"])\n",
    "\n",
    "    # Replace NaN values with the median of each column\n",
    "    X = X.fillna(X.median())\n",
    "\n",
    "    # Normalize the features\n",
    "    scaler = StandardScaler()\n",
    "    X = scaler.fit_transform(X)\n",
    "\n",
    "    # Step 2: Feature selection\n",
    "\n",
    "    # Remove features with low variance\n",
    "    selector = VarianceThreshold(threshold=0.01)\n",
    "    X = selector.fit_transform(X)\n",
    "\n",
    "    # Remove correlated features\n",
    "    def remove_correlated_features(data, threshold=0.9):\n",
    "        corr_matrix = pd.DataFrame(data).corr().abs()\n",
    "        upper = corr_matrix.where(\n",
    "            np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool_)\n",
    "        )\n",
    "        to_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n",
    "        return pd.DataFrame(data).drop(columns=to_drop, axis=1).values\n",
    "\n",
    "    X = remove_correlated_features(X)\n",
    "\n",
    "    # Step 3: Run Extremely Randomized Trees (ExtraTreesClassifier)\n",
    "    X_train, X_test, y_train, y_test = train_test_split(\n",
    "        X, y, test_size=0.2, random_state=42\n",
    "    )\n",
    "\n",
    "    clf = RandomForestClassifier(random_state=42)\n",
    "    clf.fit(X_train, y_train)\n",
    "    # Step 4: Predict probabilities and adjust threshold\n",
    "    y_pred_proba = clf.predict_proba(X_test)[:, 1]  # Probability of positive class\n",
    "\n",
    "    # Define a custom threshold to maximize accuracy for negatives\n",
    "    custom_threshold = 0.1  # Adjust this threshold as needed\n",
    "\n",
    "    y_pred = (y_pred_proba >= custom_threshold).astype(int)\n",
    "\n",
    "    # Step 5: Evaluate the model\n",
    "    report = classification_report(y_test, y_pred)\n",
    "    print(\"Classification Report:\")\n",
    "    print(report)\n",
    "\n",
    "    accuracy = accuracy_score(y_test, y_pred)\n",
    "    print(f\"Accuracy: {accuracy:.2f}\")\n",
    "\n",
    "    conf_matrix = confusion_matrix(y_test, y_pred)\n",
    "    print(\"Confusion Matrix:\")\n",
    "    print(conf_matrix)\n",
    "    # Get feature importances and indices of top features\n",
    "    importances = clf.feature_importances_\n",
    "    print(importances)\n",
    "classify(combined_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True Positives (TP): 1\n",
      "False Positives (FP): 0\n",
      "False Negatives (FN): 2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAB8YAAAGoCAYAAAA0KY9AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/aElEQVR4nO3dd5gV1d0H8N+ylKVJVxCVslgAjSQY9VURjIUIGrHRNIJColEgGktEfFViQZPYC4pJFFtexJqgRuUNJlFiRZNoRBEXNfiKWFAjSj3vHzx7w2V3YSnL7g6fz/Ps83Bn5s45d5g7Z2a+d84pSCmlAAAAAAAAAICMqlPdFQAAAAAAAACAqiQYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGIf1VFBQEBdddFF1V2Othg8fHk2aNKnuagBkTseOHWP48OG510899VQUFBTEU089tcnKqA3tDACVd/vtt0dBQUHMmzdvnctWRbuyqZ166qlx8MEHV2kZe++9d5xzzjlVWgZQu9WG4+WmsD7XBmteq9Q0zz//fNSvXz/eeeedKln/P//5z6hbt268+uqrVbJ+oObbUtqGqnDRRRdFQUFBdVejQitXroxdd901Lr300iov6+abb44ddtghlixZUuVlUT0E41SJkpKSGDVqVOy0007RqFGjaNSoUXTr1i1OO+20+Pvf/17d1atSffr0iYKCgnX+bWzosXjx4rjooos09MAWpTRcKP0rKiqKnXbaKUaNGhULFiyo7upV2qOPPir8BjKnMufANeFG1Zrn6y1btoxvf/vb8Zvf/CZWrly5Wepw0003xe23375ZytqUSkpK4le/+lWcd955EbH+1z4dO3bMm7711ltHr1694sEHH8wr56c//WnceOON8cEHH2zujwhUsTXP51f/O/fcc6u7emtVE65FZs6cGRdddFEsWrRos5S3KY0bNy6GDBkSHTp0yE1bWzsye/bsiPhP0FX6V69evejcuXOccMIJ8fbbb+fW1a1bt+jfv39ccMEFm/2zARsnC21DUVFRzJ8/v8z8Pn36xK677loNNSurNucJv/3tb+O9996LUaNG5aZVdr8pvQYZPXp0mfWWtjH33Xdfbtrw4cNj6dKlccstt1Tth6La1K3uCpA906ZNi0GDBkXdunXjuOOOi9133z3q1KkTs2fPjgceeCAmTpwYJSUleSfCWTJu3LgYOXJk7vULL7wQ1113XZx33nnRtWvX3PRvfOMbG1XO4sWLY/z48RGxqoEF2JL87Gc/i06dOsXXX38dTz/9dEycODEeffTRePXVV6NRo0abrR77779/fPXVV1G/fv31et+jjz4aN954Y7nh+FdffRV16zpFA2qfO++8M+/1HXfcEU8++WSZ6aufE1eX7bbbLiZMmBAREQsXLow77rgjRowYEW+++WZcfvnlm7Ss73//+zF48OBo0KBBbtpNN90UrVu3LvNk34a2K5vLtddeG506dYoDDjggIjbs2qdHjx5x5plnRkTE+++/H7fcckscddRRMXHixDjllFMiIuKII46IrbbaKm666ab42c9+tjk+GrCZlZ7Pr66mBAfrsjmvRda8Npg5c2aMHz8+hg8fHs2bN89b9o033og6dWrmM1CvvPJKTJ8+PWbOnFlm3upt8uq23XbbvNdjxoyJb3/727Fs2bKYNWtWTJo0KR555JH4xz/+kVv2lFNOiX79+sXcuXOjuLi4aj4MUGVqc9uwZMmSuPzyy+P666+v7qpUaG15wvnnn1+jf4Twi1/8IgYPHhzNmjUrM6+y+82tt94aY8eOLdO+rKmoqCiGDRsWV111VYwePbpGP0nPhnHXlU1q7ty5MXjw4OjQoUP87//+b7Rr1y5v/hVXXBE33XTTOk/Uv/zyy2jcuHFVVrXKrNmtYFFRUVx33XVx8MEHrzXArs2fGWBzO/TQQ2OPPfaIiIiRI0dGq1at4qqrroqHH344hgwZUmb5qjrG1qlTJ4qKijbpOjf1+gA2l+OPPz7v9bPPPhtPPvlkmelrWrx48Wb9UVNERLNmzfLqdfLJJ8fOO+8cN9xwQ1x88cVRr169TVZWYWFhFBYWVmrZqmhXNpVly5bF3XffnQuvIzbs2qd9+/Z52/6EE06ILl26xNVXX51bd506deKYY46JO+64I8aPH+9mFGTQ6ufztc36XotsjPVpE1b/AVZNc9ttt8UOO+wQe++9d5l5a7bJFenVq1ccc8wxERFx4oknxk477RRjxoyJyZMnx9ixYyMi4qCDDooWLVrE5MmT/bAKaqHa3Db06NGj0sFrTVS3bt0a+5DGyy+/HH/729/iyiuvLHd+Zfab7t27xxtvvBGXX355XHfddessc+DAgfHzn/88ZsyYEd/5znc2qN7UXDXzZ4TUWj//+c/jyy+/jNtuu61MKB6x6gA7ZsyY2H777XPTSsfDnjt3bvTr1y+aNm0axx13XESsCjLOPPPM2H777aNBgwax8847xy9/+ctIKeXeP2/evCgoKCi3K8I1uywvHSvjrbfeyv26tlmzZnHiiSfG4sWL8967ZMmSOOOMM6JNmzbRtGnT+N73vhf/+te/NnIL5dfjn//8ZwwdOjRatGgR++23X0Ss+rVWeTeRhg8fHh07dsx95jZt2kRE5G4Uldc9+/z582PAgAHRpEmTaNOmTZx11lmxYsWKTfIZAGqS0pPUkpKStbYrK1eujGuuuSa6d+8eRUVFsc0228TJJ58cn376ad76UkpxySWXxHbbbReNGjWKAw44IF577bUy5VY0ftVzzz0X/fr1ixYtWkTjxo3jG9/4Rlx77bURsep4fuONN0ZEfrfDpco7nr/88stx6KGHxlZbbRVNmjSJAw88MJ599tm8ZUq7kHrmmWfiJz/5SbRp0yYaN24cRx55ZCxcuHD9NypAFSjtSvCll16K/fffPxo1apTrlrui4YbKGzN10aJFcfrpp+euE7p06RJXXHHFBneF3qhRo9h7773jyy+/zB0z33777Tj22GOjZcuWufmPPPJImfdef/310b1792jUqFG0aNEi9thjj7jnnnty89ccY7xjx47x2muvxZ/+9KdcG1B6/r9muzJq1Kho0qRJmWuViIghQ4ZE27Zt887vH3vssejVq1c0btw4mjZtGv379y/Tfn3wwQdx4oknxnbbbRcNGjSIdu3axRFHHLHOMdCffvrp+Oijj+Kggw5a1+ZcL23bto2uXbtGSUlJ3vSDDz443nnnnXjllVc2aXlAzfbOO+/EqaeeGjvvvHM0bNgwWrVqFccee+w6j1EREXPmzImjjz462rZtG0VFRbHddtvF4MGD47PPPstb7q677oqePXtGw4YNo2XLljF48OB47733NrjOq1+LREQsX748Lr744iguLo4GDRpEx44d47zzziszVumLL74Yffv2jdatW0fDhg2jU6dOcdJJJ+Uts3rbeNFFF8XZZ58dERGdOnXKtSGrty+l7eWLL74YBQUFMXny5DL1ffzxx6OgoCCmTZuWmzZ//vw46aSTYptttokGDRpE9+7d4ze/+U2Z966rzavIQw89FN/5znc26Q+d1tzuERH16tWLPn36xMMPP7zJygGqX21oG84777xYsWJFpXufqmx5N954Y3Tu3DkaNmwYe+65Z/zlL38pkx8sXbo0LrjggujZs2c0a9YsGjduHL169YoZM2bklllXnrDmGOO77rprrpeo1a1cuTLat2+f+6FS6bTK3GurTLtXnoceeijq168f+++//zqXrUjHjh3jhBNOiFtvvTXef//9dS7fs2fPaNmypfYko2rmT0CotaZNmxZdunSJvfbaa73et3z58ujbt2/st99+8ctf/jIaNWoUKaX43ve+FzNmzIgRI0ZEjx494vHHH4+zzz475s+fH1dfffUG13PgwIHRqVOnmDBhQsyaNSt+9atfxdZbbx1XXHFFbpmRI0fGXXfdFUOHDo199tkn/vjHP0b//v03uMzyHHvssbHjjjvGZZddlhf2r0ubNm1i4sSJ8aMf/SiOPPLIOOqooyIiv4vCFStWRN++fWOvvfaKX/7ylzF9+vS48soro7i4OH70ox9t0s8BUN3mzp0bERGtWrWKiPLblYhVTwTefvvtceKJJ8aYMWOipKQkbrjhhnj55ZfjmWeeyT0heMEFF8Qll1wS/fr1i379+sWsWbPikEMOiaVLl66zLk8++WQcdthh0a5du/jxj38cbdu2jddffz2mTZsWP/7xj+Pkk0+O999/v9zuhcvz2muvRa9evWKrrbaKc845J+rVqxe33HJL9OnTJ/70pz+VaXNHjx4dLVq0iAsvvDDmzZsX11xzTYwaNSqmTJmyXtsUoKp8/PHHceihh8bgwYPj+OOPj2222Wa93r948eLo3bt3zJ8/P04++eTYYYcdYubMmTF27Nj4v//7v7jmmms2qF5vv/12FBYWRvPmzWPBggWxzz77xOLFi2PMmDHRqlWrmDx5cnzve9+L++67L4488siIWNUd35gxY+KYY46JH//4x/H111/H3//+93juuedi6NCh5ZZzzTXXxOjRo6NJkyYxbty4iIgKt8GgQYPixhtvjEceeSSOPfbYvG3w+9//PoYPH557Gv3OO++MYcOGRd++feOKK66IxYsXx8SJE2O//faLl19+Ofcj26OPPjpee+21GD16dHTs2DE+/PDDePLJJ+Pdd9/NLVOemTNnRkFBQXzzm99c3027VsuWLYv33nsv14aX6tmzZ0REPPPMM5u8TKD6ffbZZ/HRRx/lTWvdunW88MILMXPmzBg8eHBst912MW/evJg4cWL06dMn/vnPf1bYw8jSpUujb9++sWTJkhg9enS0bds25s+fH9OmTYtFixblul699NJL47//+79j4MCBMXLkyFi4cGFcf/31sf/++8fLL79cpnvyyljzWmTkyJExefLkOOaYY+LMM8+M5557LiZMmBCvv/56PPjggxER8eGHH8YhhxwSbdq0iXPPPTeaN28e8+bNiwceeKDCco466qh4880347e//W1cffXV0bp164iIXNCxuj322CM6d+4c9957bwwbNixv3pQpU6JFixbRt2/fiIhYsGBB7L333lFQUBCjRo2KNm3axGOPPRYjRoyIzz//PE4//fSI2LA2L2JV6P7uu+/Gt771rXLnr1ixosy+UFRUFE2aNKlwnRFlt3upnj17xsMPPxyff/55bLXVVmtdB1Cz1Oa2oVOnTrng9dxzz13rU+OVLW/ixIkxatSo6NWrV5xxxhkxb968GDBgQLRo0SK222673Po+//zz+NWvfhVDhgyJH/zgB/HFF1/Er3/96+jbt288//zz0aNHj0rlCasbNGhQXHTRRfHBBx9E27Ztc9OffvrpeP/992Pw4MG5aZW517Yh7V6pmTNnxq677lphr14V7TdrGjduXNxxxx2Vfmr8W9/6VjzzzDPrXI5aKMEm8tlnn6WISAMGDCgz79NPP00LFy7M/S1evDg3b9iwYSki0rnnnpv3noceeihFRLrkkkvyph9zzDGpoKAgvfXWWymllEpKSlJEpNtuu61MuRGRLrzwwtzrCy+8MEVEOumkk/KWO/LII1OrVq1yr1955ZUUEenUU0/NW27o0KFl1rkuU6dOTRGRZsyYUaYeQ4YMKbN87969U+/evctMHzZsWOrQoUPu9cKFCyusS+k2/dnPfpY3/Zvf/Gbq2bNnpesOUNPcdtttKSLS9OnT08KFC9N7772X/ud//ie1atUqNWzYMP3rX/+qsF35y1/+kiIi3X333XnT//CHP+RN//DDD1P9+vVT//7908qVK3PLnXfeeSki0rBhw3LTZsyYkXeMX758eerUqVPq0KFD+vTTT/PKWX1dp512WqroNGzNY/uAAQNS/fr109y5c3PT3n///dS0adO0//77l9k2Bx10UF5ZZ5xxRiosLEyLFi0qtzyAqlLesa53794pItLNN99cZvmKzm07dOiQd+y9+OKLU+PGjdObb76Zt9y5556bCgsL07vvvrvWevXu3TvtsssuuWuT119/PY0ZMyZFRDr88MNTSimdfvrpKSLSX/7yl9z7vvjii9SpU6fUsWPHtGLFipRSSkcccUTq3r37WssrPT6XlJTkpnXv3r3cc/4125WVK1em9u3bp6OPPjpvuXvvvTdFRPrzn/+cq1vz5s3TD37wg7zlPvjgg9SsWbPc9E8//TRFRPrFL36x1jqX5/jjj8+7ZipPedc+q+vQoUM65JBDctv+b3/7Wxo8eHCKiDR69Ogyy9evXz/96Ec/Wu+6AjVX6TGxvL+UUt79olJ//etfU0SkO+64IzdtzePlyy+/nCIiTZ06tcKy582blwoLC9Oll16aN/0f//hHqlu3bpnpFdV9bdcipfeTRo4cmffes846K0VE+uMf/5hSSunBBx9MEZFeeOGFtZa5Ztv4i1/8okybUmrN9nLs2LGpXr166ZNPPslNW7JkSWrevHnefbERI0akdu3apY8++ihvfYMHD07NmjXL/Z9Ups0rz/Tp01NEpN///vdl5pWeF6z5V941129+85u0cOHC9P7776dHHnkkdezYMRUUFJTZhvfcc0+KiPTcc8+td12B6pGFtuGFF15Ic+fOTXXr1k1jxozJze/du3fesbOy5S1ZsiS1atUqffvb307Lli3LLXf77beniMi7lli+fHlasmRJ3vo+/fTTtM022+Qd79eWJ5TmFaXeeOONFBHp+uuvz1vu1FNPTU2aNMn9n1T2Xltl273ybLfddmWuh1Ja935TqkOHDql///4ppZROPPHEVFRUlN5///2U0n/2mfL2kR/+8IepYcOG611faj5dqbPJfP755xER5f6is0+fPtGmTZvcX2kXsqtb8ynmRx99NAoLC2PMmDF5088888xIKcVjjz22wXVdfVy8iFXjFH388ce5z/Doo49GRJQpu/RXspvKmvXY1Mr7nG+//XaVlgmwORx00EHRpk2b2H777WPw4MHRpEmTePDBB6N9+/a5ZdZsV6ZOnRrNmjWLgw8+OD766KPcX8+ePaNJkya5LqamT58eS5cujdGjR+d1I1WZNuDll1+OkpKSOP3008v8onhDug1csWJFPPHEEzFgwIDo3Llzbnq7du1i6NCh8fTTT+farlI//OEP88rq1atXrFixIt555531Lh+gKjRo0CBOPPHEDX7/1KlTo1evXtGiRYu84/lBBx0UK1asiD//+c/rXMfs2bNz1yZdu3aN66+/Pvr375/rNvbRRx+NPffcMzfcUcSq65wf/vCHMW/evPjnP/8ZERHNmzePf/3rX/HCCy9s8OdZm4KCgjj22GPj0UcfjX//+9+56VOmTIn27dvn6vfkk0/GokWLYsiQIXnbpLCwMPbaa69cG9ewYcOoX79+PPXUU2W6NlyXjz/+OFq0aLHRn+mJJ57Ibfvdd989pk6dGt///vfzeu8qVfp/DGTPjTfeGE8++WTeX8Sq41SpZcuWxccffxxdunSJ5s2bx6xZsypcX+lTf48//ni5w09ERDzwwAOxcuXKGDhwYN6xsm3btrHjjjvmdTm7Nmu7Fim9n/STn/wk7z1nnnlmRERuSI7Sa4Vp06bFsmXLKlXu+ho0aFAsW7Ys72m8J554IhYtWhSDBg2KiFVDSN1///1x+OGHR0opb7v07ds3Pvvss9x239A27+OPP46IqLAN6dixY5l94Zxzzimz3EknnRRt2rSJbbfdNvr37x9ffvllTJ48ucy4sqXlaD+g9qnNbUNEROfOneP73/9+TJo0Kf7v//5vo8p78cUX4+OPP44f/OAHeWN/H3fccWWOp4WFhVG/fv2IWNWt+SeffBLLly+PPfbYY63bZ2122mmn6NGjR17vgytWrIj77rsvDj/88Nz/SWXvtW1Mu7eu65CK9pvynH/++bF8+fJKdXnfokWL+Oqrryrcd6i9dKXOJtO0adOIiLwbNqVuueWW+OKLL2LBggVx/PHHl5lft27dvO4/IlaNHbLtttvm1luqa9euufkbaocddsh7XXpg/fTTT2OrrbaKd955J+rUqRPFxcV5y+28884bXGZ5OnXqtEnXt7qioqIy3Wm1aNFivW+AAdREN954Y+y0005Rt27d2GabbWLnnXeOOnX+83u/8tqVOXPmxGeffRZbb711uev88MMPI+I/7cuOO+6YN79NmzbrDARKu/Pbdddd1+8DVWDhwoWxePHictufrl27xsqVK+O9996L7t2756avrY0DqAnat2+fu3GzIebMmRN///vfy+06NuI/x/O16dixY9x6661RUFAQRUVFseOOO+a1D++88065w0Otfi2y6667xk9/+tOYPn167LnnntGlS5c45JBDYujQobHvvvtu4Kcra9CgQXHNNdfE7373uxg6dGj8+9//jkcffTROPvnk3A+h5syZExH/GW91TaVdyTZo0CCuuOKKOPPMM2ObbbaJvffeOw477LA44YQT8rpIrEhaj+GfKrLXXnvFJZdcEgUFBdGoUaPo2rVrhd1TppQ26Xi0QM2x5557lgk0IyK++uqrmDBhQtx2220xf/78vOPOmuPBrq5Tp07xk5/8JK666qq4++67o1evXvG9730vjj/++FwwMmfOnEgplTnPL1VRF61rWtu1SOn9pC5duuS9p23bttG8efPctUbv3r3j6KOPjvHjx8fVV18dffr0iQEDBsTQoUOjQYMGlarHuuy+++6xyy67xJQpU2LEiBERseqHVa1bt861FwsXLoxFixbFpEmTYtKkSeWup7Rd3dg2r6I2pHHjxnHQQQet8/0XXHBB9OrVKwoLC6N169bRtWvXvLBozXK0H1D71Oa2odT5558fd955Z1x++eVx7bXXlplf2fJK24s125O6deuWO/zR5MmT48orr4zZs2fnBc8bkz8MGjQozjvvvJg/f360b98+nnrqqfjwww9zP64q/TyVude2se3e2q5DKtpvyrP6jxfOPffcSpWpPckewTibTLNmzaJdu3bx6quvlplXelNp3rx55b63QYMGeYHG+qjowLRixYoK31M6Dt+aNsWNnvWx+q/dShUUFJRbj7V9nvJU9BkBsmBdJ73ltSsrV66MrbfeOu6+++5y31NRwFLb1JQ2DqAi5Z0Dr82a58ErV66Mgw8+uNynySJWPd2wLpW9Cb8uXbt2jTfeeCOmTZsWf/jDH+L++++Pm266KS644IIYP378Rq8/ImLvvfeOjh07xr333htDhw6N3//+9/HVV1/l3ZBauXJlRKwaZ7y8gHv14OD000+Pww8/PB566KF4/PHH47//+79jwoQJ8cc//nGtY3m3atVqk/zIqnXr1pXe9osWLSp3fEAgu0aPHh233XZbnH766fFf//Vf0axZsygoKIjBgwfnjnUVufLKK2P48OHx8MMPxxNPPBFjxoyJCRMmxLPPPhvbbbddrFy5MgoKCuKxxx4r95x5XWNal6rMDfh13UQvKCiI++67L5599tn4/e9/H48//nicdNJJceWVV8azzz5b6bqsy6BBg+LSSy+Njz76KJo2bRq/+93vYsiQIbl2oXSbHn/88WXGIi9VOv7shrZ5pWOAb2wbsttuu1Wq/SgtR/sB2VEb2oZSnTt3juOPP77C4HVTlxcRcdddd8Xw4cNjwIABcfbZZ8fWW28dhYWFMWHChNwDHBti0KBBMXbs2Jg6dWqcfvrpce+990azZs3iu9/9bt7nqcy9to1p9zbVdUipcePGxZ133hlXXHFFDBgwoMLlPv3002jUqNF6X79S8wnG2aT69+8fv/rVr+L555+PPffcc6PW1aFDh5g+fXp88cUXeU+Nz549Ozc/4j9Pwi1atCjv/RvzRHmHDh1i5cqVMXfu3Lyn9N54440NXmdltWjRotzuztf8PH6pBLB+iouLY/r06bHvvvuu9aS2tH2ZM2dOXvflCxcuXOeJeGlPI6+++upab9pU9hjepk2baNSoUbntz+zZs6NOnTqx/fbbV2pdADVdixYtypzTL126tEw3hMXFxfHvf/97kwTbFenQoUOFx97S+aUaN24cgwYNikGDBsXSpUvjqKOOiksvvTTGjh0bRUVF5a5/fc/lBw4cGNdee218/vnnMWXKlOjYsWPsvffeufml7c/WW29dqe1SXFwcZ555Zpx55pkxZ86c6NGjR1x55ZVx1113VfieXXbZJe6+++747LPPck/YVKX58+fH0qVLc0/pA1uG++67L4YNGxZXXnllbtrXX39dpn2oyG677Ra77bZbnH/++TFz5szYd9994+abb45LLrkkiouLI6UUnTp1qtSPqDZE6f2kOXPm5B2/FixYEIsWLcprPyJW/fhp7733jksvvTTuueeeOO644+J//ud/YuTIkeWuf33bj0GDBsX48ePj/vvvj2222SY+//zzGDx4cG5+mzZtomnTprFixYpKtR8b0ubtsssuERFRUlKyXnXfUCUlJVGnTp0q+z8GNr/a1jacf/75cdddd5U7VFBlyyttL95666044IADctOXL18e8+bNy/1oKWLV9uncuXM88MADee3EhRdemLfO9W1DOnXqFHvuuWdMmTIlRo0aFQ888EAMGDAg7wnvyt5rK7W+7V7EqnZkU7YhxcXFcfzxx8ctt9xSbi9hpUpKSlyLZJQxxtmkzjnnnGjUqFGcdNJJsWDBgjLz1+dptX79+sWKFSvihhtuyJt+9dVXR0FBQRx66KERsapbwNatW5cZS/Cmm27agE+wSum6r7vuurzp11xzzQavs7KKi4tj9uzZsXDhwty0v/3tb/HMM8/kLdeoUaOIKPuDAADKN3DgwFixYkVcfPHFZeYtX748dzw96KCDol69enH99dfntVuVaQO+9a1vRadOneKaa64pc3xefV2NGzeOiHUfwwsLC+OQQw6Jhx9+OK/XlQULFsQ999wT++23X657XIDarri4uMw5/aRJk8o8MT5w4MD461//Go8//niZdSxatCiWL1++0XXp169fPP/88/HXv/41N+3LL7+MSZMmRceOHaNbt24R8Z9xU0vVr18/unXrFimltY6d17hx4/U6jx80aFAsWbIkJk+eHH/4wx9i4MCBefP79u0bW221VVx22WXlllt6bbF48eL4+uuv8+YVFxdH06ZNY8mSJWutw3/9139FSileeumlStd7Y5SWs88++2yW8oCaobCwsMy9o+uvv36dveh9/vnnZY7/u+22W9SpUyd3fDvqqKOisLAwxo8fX6aMlFKZY/qG6NevX0SUvXa46qqrImLVAyURq55CW7MOPXr0iIhY6/G4stcRpbp27Rq77bZbTJkyJaZMmRLt2rWL/fffPze/sLAwjj766Lj//vvL7QFy9XtTG9rmtW/fPrbffvt48cUXK1XnjfXSSy9F9+7dN8uPuIDNo7a1DasHrx988EHevMqWt8cee0SrVq3i1ltvzfsMd999d5mHNkqfPF99fc8991zetUzEhuUJgwYNimeffTZ+85vfxEcffZTXa1VE5e+1bWi7F7HqOuTVV19d53Lr4/zzz49ly5bFz3/+8wqXmTVrlmuRjPLEOJvUjjvuGPfcc08MGTIkdt555zjuuONi9913j5RSlJSUxD333BN16tQpM+5reQ4//PA44IADYty4cTFv3rzYfffd44knnoiHH344Tj/99Lzxv0eOHBmXX355jBw5MvbYY4/485//HG+++eYGf44ePXrEkCFD4qabborPPvss9tlnn/jf//3feOuttzZ4nZV10kknxVVXXRV9+/aNESNGxIcffhg333xzdO/ePT7//PPccg0bNoxu3brFlClTYqeddoqWLVvGrrvuusnGtQXImt69e8fJJ58cEyZMiFdeeSUOOeSQqFevXsyZMyemTp0a1157bRxzzDHRpk2bOOuss2LChAlx2GGHRb9+/eLll1+Oxx57bJ3d8dWpUycmTpwYhx9+ePTo0SNOPPHEaNeuXcyePTtee+21XIjTs2fPiIgYM2ZM9O3bNwoLC/Oe3FjdJZdcEk8++WTst99+ceqpp0bdunXjlltuiSVLlqz1BB6gthk5cmSccsopcfTRR8fBBx8cf/vb3+Lxxx8vc+w9++yz43e/+10cdthhMXz48OjZs2d8+eWX8Y9//CPuu+++mDdv3kZ3n3ruuefGb3/72zj00ENjzJgx0bJly5g8eXKUlJTE/fffnxuu45BDDom2bdvGvvvuG9tss028/vrrccMNN0T//v3zer1aU8+ePWPixIlxySWXRJcuXWLrrbeucHzwiFU/vOrSpUuMGzculixZUuaG1FZbbRUTJ06M73//+/Gtb30rBg8eHG3atIl33303Hnnkkdh3333jhhtuiDfffDMOPPDAGDhwYHTr1i3q1q0bDz74YCxYsKDCdqjUfvvtF61atYrp06evta6bypNPPhk77LDDWrt3B7LnsMMOizvvvDOaNWsW3bp1i7/+9a8xffr0XHfcFfnjH/8Yo0aNimOPPTZ22mmnWL58edx555254DdiVVBxySWXxNixY2PevHkxYMCAaNq0aZSUlMSDDz4YP/zhD+Oss87aqPrvvvvuMWzYsJg0aVIsWrQoevfuHc8//3xMnjw5BgwYkHvqb/LkyXHTTTfFkUceGcXFxfHFF1/ErbfeGltttVUuXC9P6XXEuHHjYvDgwVGvXr04/PDDc4F5eQYNGhQXXHBBFBUVxYgRI8oMOXX55ZfHjBkzYq+99oof/OAH0a1bt/jkk09i1qxZMX369Pjkk08iYsPbvIiII444Ih588MFIKVVpD4jLli2LP/3pT3HqqadWWRnA5lcb24bS7rrfeOON6N69e256ZcurX79+XHTRRTF69Oj4zne+EwMHDox58+bF7bffHsXFxXnH0sMOOyweeOCBOPLII6N///5RUlISN998c3Tr1i3+/e9/55bbkDxh4MCBcdZZZ8VZZ50VLVu2LNO7SGXvtW1ouxexqg25+OKL409/+lMccsgh6/tfUa7SHy9Mnjy53PkvvfRSfPLJJ3HEEUdskvKoYRJUgbfeeiv96Ec/Sl26dElFRUWpYcOGaZdddkmnnHJKeuWVV/KWHTZsWGrcuHG56/niiy/SGWeckbbddttUr169tOOOO6Zf/OIXaeXKlXnLLV68OI0YMSI1a9YsNW3aNA0cODB9+OGHKSLShRdemFvuwgsvTBGRFi5cmPf+2267LUVEKikpyU376quv0pgxY1KrVq1S48aN0+GHH57ee++9Mutcl6lTp6aISDNmzFhnPUrdddddqXPnzql+/fqpR48e6fHHH0/Dhg1LHTp0yFtu5syZqWfPnql+/fp59apom5aWC1BblR6vX3jhhQqXWVu7klJKkyZNSj179kwNGzZMTZs2Tbvttls655xz0vvvv59bZsWKFWn8+PGpXbt2qWHDhqlPnz7p1VdfTR06dEjDhg3LLTdjxowyx/iUUnr66afTwQcfnJo2bZoaN26cvvGNb6Trr78+N3/58uVp9OjRqU2bNqmgoCDv2FxeOzNr1qzUt2/f1KRJk9SoUaN0wAEHpJkzZ1Zq21RUR4Cqdtppp5U59+zdu3fq3r17ucuvWLEi/fSnP02tW7dOjRo1Sn379k1vvfVWmWNvSquuE8aOHZu6dOmS6tevn1q3bp322Wef9Mtf/jItXbp0rfVaWx1WN3fu3HTMMcek5s2bp6KiorTnnnumadOm5S1zyy23pP333z+1atUqNWjQIBUXF6ezzz47ffbZZ7llyrvW+OCDD1L//v1T06ZNU0Sk3r17p5TWfsweN25ciojUpUuXCus8Y8aM1Ldv39SsWbNUVFSUiouL0/Dhw9OLL76YUkrpo48+SqeddlraZZddUuPGjVOzZs3SXnvtle699951bo+UUhozZsxayy/v2md1HTp0SP37919nOStWrEjt2rVL559/fqXqBdQe6zqf//TTT9OJJ56YWrdunZo0aZL69u2bZs+evc7z8LfffjuddNJJqbi4OBUVFaWWLVumAw44IE2fPr1MGffff3/ab7/9UuPGjVPjxo3TLrvskk477bT0xhtvbFTdSy1btiyNHz8+derUKdWrVy9tv/32aezYsenrr7/OLTNr1qw0ZMiQtMMOO6QGDRqkrbfeOh122GG543Wp8q4NLr744tS+fftUp06dvPalvPYypZTmzJmTIiJFRHr66afLrfOCBQvSaaedlrbffvtUr1691LZt23TggQemSZMm5ZapTJtXkVmzZqWISH/5y1/yplemTS79v546deo6y3nsscdSRKQ5c+asc1mg5shq2zBs2LAUEeUe5ypb3nXXXZc6dOiQGjRokPbcc8/0zDPPpJ49e6bvfve7uWVWrlyZLrvsstxy3/zmN9O0adPWK09YW26w7777pohII0eOrHA7rOteW2XbvYp84xvfSCNGjMibVtl2uaJrkDlz5qTCwsJy25if/vSnaYcddiiTQ5ENBSmtR9/WAAAAwBbp7bffjl122SUee+yxOPDAA6usnIceeiiGDh0ac+fOjXbt2lVZOQBsPgceeGBsu+22ceedd1ZZGQMGDIiCgoJ48MEHq6wMgOq0cuXKaNOmTRx11FFx6623Vnd1Nps777wzTjvttHj33XejefPmVVrWkiVLomPHjnHuuefGj3/84yoti+phjHEAAABgnTp37hwjRoyIyy+/vErLueKKK2LUqFFCcYAMueyyy2LKlCnxzjvvVMn6X3/99Zg2bVq549wC1EZff/11mXG577jjjvjkk0+iT58+1VOpanLcccfFDjvsEDfeeGOVl3XbbbdFvXr14pRTTqnysqgenhgHAAAAAACAGuKpp56KM844I4499tho1apVzJo1K379619H165d46WXXor69etXdxWhVqpb3RUAAAAAAAAAVunYsWNsv/32cd1118Unn3wSLVu2jBNOOCEuv/xyoThsBE+MAwAAAAAAAJBpxhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEyrW9kFCwoKqrIeANRAKaUNep82A2DLsqHtRYQ2A2BLo80AoLK0GQBUVmXbDE+MAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKbVre4KAAA1V0op9++CgoJqrAkAAAAAAGw4T4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYZYxwAMmz1McKre13GKAcAAAAAoLp4YhwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMs0Y41ugNceINeYrm5p9DDavTTmOeFVaWz0dJwBqv9WP847rAAAA1HSyjC2PJ8YBAAAAAAAAyDTBOAAAAAAAAACZpiv1WmpTdpu7MevSrUR22ceg5qotXaevD90WAQAAAADrS5bB+hCMAwAAUGNuJkS4oQAAAABserpSBwAAAAAAACDTBOMAAAAAAAAAZJqu1Guo2jJ+7LrqqQvEmss+BrVHbfm+bkrGHAcAAAAAasu9UVlG7SAYBwAA2EJk4YaCmwkAAADAhtCVOgAAAAAAAACZJhgHAAAAAAAAINN0pV5D1JYuDdeXMWJrDvsY1B5Z/b5uDN91AAAAANgyZPH+qPubNYNgHAAAIKPcTAAAAABYRVfqAAAAAAAAAGSaYBwAAAAAAACATNOVejXJYpeGlbH659blYdWyj61iP6M22FK/rxvDdx0AAAAAsmFLvD/q/mb1EIwDAABkhJsJbiYAAAAA5dOVOgAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADKtbnVXYEuRUqruKtQ4a26TgoKCaqpJdtjPylp9m9jHqEl8Xzct33VgS6U9Kct1BgAAQM3mWrYs17KbhyfGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpdau7AlmWUqruKtQqa26vgoKCaqpJ7WEfWz/2MaqT7+vm47sOAAAAADWHe6Prb/Vt5v7mpiMYBwAAqCXcTFh/fjAFAAAAROhKHQAAAAAAAICME4wDAAAAAAAAkGm6Uq9Cq3fRp8vDddOl4fpbc5vZz9bOPkZ18n3dfHzXgSzTnqw/7QIAAED1ci27/lzLVg1PjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0+pWdwW2FAUFBXmvU0rVVJOaY81twsZbfZvax1axn1FT+b5uWr7rAAAAAFA7yMzKcn9z8xCMAwAA1FJuJpTlZgIAAABQHl2pAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINPqVncFtlQFBQV5r1NK1VSTzWvNz03VsY9B7bGlfl83hu86QPm2xDZFmwAAAFC7uZZlc/HEOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmjHGa4isjp9gjISawz4GtUdWv68bw3cdYMNksU3RJgAAAGTb6td9WbiOjXAtW1N4YhwAAAAAAACATBOMAwAAAAAAAJBpulKvodbVpUJN6TpC1w+1l30Mao8sdoO7Lr77AAAAAIAsg01JMA4AALCFWNuFek25mRDhhgIAAACw6elKHQAAAAAAAIBME4wDAAAAAAAAkGm6Uq+lNqZrwTW7SNRNIeWxj0HNlcUxxx0nAAAAAID1JctgfQjGAQAA2OgbAKvfUHAzAQAAAKhpdKUOAAAAAAAAQKYJxgEAAAAAAADINF2pb4F0a0hVs4/B5rW271xNGn/csQEAAAAAqCncr9zyCMYBAADYaG4oAAAAADWZrtQBAAAAAAAAyDTBOAAAAAAAAACZpit1AMiwje3WdvUxynWRCwAAAABAbeWJcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINGOMAwAVMq44AAAAAABZ4IlxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpdau7AlDTpZTyXhcUFFRTTQAAAAAAAMonz4C188Q4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaMcbZIqw5rkZ1rsuYHgAAAAAAQHlqSp4hyyCLPDEOAAAAAAAAQKYJxgEAAAAAAADINF2pkxmbsnuRqrS2euqaBAAAAAAAsisLWUaEPIPayRPjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpxhin1qot43CsjzU/kzE6AAAAAACg9spilhEhz6B28sQ4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaMcapNbI6DsfaGKMDAAAAAABqjy0xy4jI/9yyDGoqT4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGRa3equAFQkpVTdVahx1twmBQUF1VQTAAAAAAAgQp6xJlkGNZUnxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaXWruwJQKqVU3VWoddbcZgUFBdVUEwAAAAAA2DLIM9aPLIOawhPjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpxhinxlhzTAljdKybcTgAAAAAAGDzkmesH1kGNYUnxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaXWruwJQkYKCgrzXKaVqqknNseY2AQAAAAAAqtfq9+5lGbIMai5PjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZFrd6q4AVFZBQUHe65RSNdVk81nzMwMAAAAAADXXlphlRMgzqB08MQ4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmWaMcWqtLI7TYQwOAAAAAADIjixmGRHyDGonT4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYZY5zMWNt4FjVpzA7jbgAAAAAAwJZpXRlBTckzZBlkkSfGAQAAAAAAAMg0wTgAAAAAAAAAmaYrdbYIG9Plx5rdlug+BAAAAAAAqAryDKg6nhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwzxjisgzE4AAAAAACAmk6eAWvniXEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpgnEAAAAAAAAAMk0wDgAAAAAAAECmCcYBAAAAAAAAyDTBOAAAAAAAAACZJhgHAAAAAAAAINME4wAAAAAAAABkmmAcAAAAAAAAgEwTjAMAAAAAAACQaYJxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAyTTAOAAAAAAAAQKYJxgEAAAAAAADINME4AAAAAAAAAJkmGAcAAAAAAAAg0wTjAAAAAAAAAGSaYBwAAAAAAACATBOMAwAAAAAAAJBpBSmlVN2VAAAAAAAAAICq4olxAAAAAAAAADJNMA4AAAAAAABApgnGAQAAAAAAAMg0wTgAAAAAAAAAmSYYBwAAAAAAACDTBOMAAAAAAAAAZJpgHAAAAAAAAIBME4wDAAAAAAAAkGmCcQAAAAAAAAAy7f8BLg2ZefZVoVAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 2000x500 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from scipy.ndimage import label, find_objects\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Example ground truth and prediction arrays\n",
    "gt = np.zeros((100, 100), dtype=np.uint8)\n",
    "pred = np.zeros((100, 100), dtype=np.uint8)\n",
    "\n",
    "# Create three blobs in the ground truth\n",
    "cv2.circle(gt, (30, 30), 10, 255, -1)\n",
    "cv2.circle(gt, (70, 30), 10, 255, -1)\n",
    "cv2.circle(gt, (50, 70), 10, 255, -1)\n",
    "\n",
    "# Create one blob in the prediction\n",
    "cv2.circle(pred, (30, 30), 12, 255, -1)  # Slightly larger or smaller than GT\n",
    "\n",
    "# Label the blobs\n",
    "gt_labels, num_gt = label(gt)\n",
    "pred_labels, num_pred = label(pred)\n",
    "\n",
    "# Find bounding boxes of blobs\n",
    "gt_slices = find_objects(gt_labels)\n",
    "pred_slices = find_objects(pred_labels)\n",
    "\n",
    "# Create empty masks for TP, FP, FN\n",
    "tp_mask = np.zeros_like(gt)\n",
    "fp_mask = np.zeros_like(gt)\n",
    "fn_mask = np.zeros_like(gt)\n",
    "\n",
    "# Calculate TP, FP, FN\n",
    "tp, fp, fn = 0, 0, 0\n",
    "\n",
    "\n",
    "# Helper function to determine if two blobs overlap\n",
    "def overlap(slice1, slice2):\n",
    "    return not (\n",
    "        slice1[0].stop < slice2[0].start\n",
    "        or slice1[0].start > slice2[0].stop\n",
    "        or slice1[1].stop < slice2[1].start\n",
    "        or slice1[1].start > slice2[1].stop\n",
    "    )\n",
    "\n",
    "\n",
    "matched_gt = set()\n",
    "matched_pred = set()\n",
    "\n",
    "for i, gt_slice in enumerate(gt_slices):\n",
    "    gt_blob = gt_labels == (i + 1)\n",
    "    match_found = False\n",
    "    for j, pred_slice in enumerate(pred_slices):\n",
    "        if overlap(gt_slice, pred_slice):\n",
    "            pred_blob = pred_labels == (j + 1)\n",
    "            tp_mask[gt_blob & pred_blob] = 255  # True positive overlap region\n",
    "            tp += 1\n",
    "            matched_gt.add(i)\n",
    "            matched_pred.add(j)\n",
    "            match_found = True\n",
    "            break\n",
    "    if not match_found:\n",
    "        fn_mask[gt_blob] = 255\n",
    "\n",
    "for j, pred_slice in enumerate(pred_slices):\n",
    "    if j not in matched_pred:\n",
    "        pred_blob = pred_labels == (j + 1)\n",
    "        fp_mask[pred_blob] = 255\n",
    "\n",
    "fn = num_gt - tp\n",
    "fp = num_pred - tp\n",
    "\n",
    "print(f\"True Positives (TP): {tp}\")\n",
    "print(f\"False Positives (FP): {fp}\")\n",
    "print(f\"False Negatives (FN): {fn}\")\n",
    "\n",
    "# Plot the images\n",
    "fig, axes = plt.subplots(1, 5, figsize=(20, 5))\n",
    "\n",
    "axes[0].imshow(gt, cmap=\"gray\")\n",
    "axes[0].set_title(\"Ground Truth\")\n",
    "\n",
    "axes[1].imshow(pred, cmap=\"gray\")\n",
    "axes[1].set_title(\"Prediction\")\n",
    "\n",
    "axes[2].imshow(tp_mask, cmap=\"gray\")\n",
    "axes[2].set_title(\"True Positives (TP)\")\n",
    "\n",
    "axes[3].imshow(fp_mask, cmap=\"gray\")\n",
    "axes[3].set_title(\"False Positives (FP)\")\n",
    "\n",
    "axes[4].imshow(fn_mask, cmap=\"gray\")\n",
    "axes[4].set_title(\"False Negatives (FN)\")\n",
    "\n",
    "for ax in axes:\n",
    "    ax.axis(\"off\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "aimiv2",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}