exploration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ranking Brain's Zones for Dementia Detection\n",
    "\n",
    "This iPython Notebook is used to create a simple ranking of the various brain's zones for dementia detection using our neural network. First, we must install the tools for dealing with NiFTY files and all the additional required libraries. We use Pipenv to deal with virtual environments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32mCourtesy Notice\u001b[0m: Pipenv found itself running within a virtual environment, so it will automatically use that environment, instead of creating its own for any project. You can set \u001b[1mPIPENV_IGNORE_VIRTUALENVS=1\u001b[0m to force pipenv to ignore that environment and create its own instead. You can set \u001b[1mPIPENV_VERBOSITY=-1\u001b[0m to suppress this warning.\n",
      "\u001b[1mInstalling dependencies from Pipfile.lock (ab9b16)...\u001b[0m\n",
      "  🐍   \u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m\u001b[1m▉\u001b[0m 0/0 — \u001b[30m00:00:00\u001b[0m\n",
      "To activate this project's virtualenv, run \u001b[33mpipenv shell\u001b[0m.\n",
      "Alternatively, run a command inside the virtualenv with \u001b[33mpipenv run\u001b[0m.\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!pipenv install"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nibabel as nib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploration of the Atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas = nib.load(\"atlas.nii.gz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get a frame from the atlas (converted to numpy array for easier access)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x12cb4f280>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAGgCAYAAABlriQwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk/klEQVR4nO3df3RU9Z3/8dfk1yQKmUhcElKJREoPVGlFfsSIZ9fWnKUtWiycbu1iG39sXWuoAnuKUoU9ohi1eyq164/KWsQtaMs5KlVbPT2x6LLlZwRXlgq0sMLRJrRqZlAkhOTz/cN+p84kznAzP+57Zp6Pc+Yc587N5JOP4otPXvdzJ+CccwIAACYV+T0AAADw8QhqAAAMI6gBADCMoAYAwDCCGgAAwwhqAAAMI6gBADCMoAYAwDCCGgAAwwhqAAAMy1hQ33///RozZozKy8vV2NiorVu3ZupbAQCQtwKZuNf3z372M33zm9/UQw89pMbGRq1YsULr1q3Tnj17NHLkyIRf29/fr7feekvDhw9XIBBI99AAAPCdc05HjhxRXV2dioqSrJldBkybNs21trZGn/f19bm6ujrX1taW9GsPHTrkJPHgwYMHDx55/zh06FDSXCxRmh0/flwdHR1avHhx9FhRUZGam5u1adOmAef39PSop6cn+tz9ZYF/ob6kEpWme3gAAPjuhHq1Ub/U8OHDk56b9qD+85//rL6+PtXU1MQcr6mp0euvvz7g/La2Nt12222DDKxUJQGCGgCQhz5ck55Uxev7Vd+LFy9WOByOPg4dOuT3kAAAMCPtK+rTTz9dxcXF6urqijne1dWl2traAecHg0EFg8F0DwMAgLyQ9hV1WVmZJk+erPb29uix/v5+tbe3q6mpKd3fDgCAvJb2FbUkLVy4UC0tLZoyZYqmTZumFStW6P3339dVV12ViW8HAEDeykhQf+1rX9Of/vQnLV26VJ2dnTr33HP1/PPPD7jADAAAJJaRG56kIhKJKBQK6SLN4qpvAEBeOuF6tUHrFQ6HVVlZmfBc36/6BgAAH4+gBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDCGoAAAwjqAEAMIygBgDAMIIaAADDSvweAJBPlu5/Jeb5MVfq00j8c8/YiX4PAcgrrKgBADCMoAYAwDCCGgAAw+iogb+I75eHohA76XiL/vBazPNel/h/M/d+ckImhwPkPFbUAAAYRlADAGAYQQ0AgGF01ChY6eikkVxp4ETM8/jOesHvf5fw6+mwUehYUQMAYBhBDQCAYQQ1AACG0VEjb9FB25Sss46XrMOOR6eNfMOKGgAAwwhqAAAMI6gBADCMjho5iw4ag6HTRr5hRQ0AgGEENQAAhhHUAAAYRkcNs+ig81OyfdPZ5rXTlui1kV2sqAEAMIygBgDAMIIaAADDbJVFKCh00IXBWiedDnyGNrKJFTUAAIYR1AAAGEZQAwBgWP6VRzCLThoAvGNFDQCAYQQ1AACGEdQAABhGR42MKYRO+pgr9XsI5uTjvul0Ozq7MaWvP+XJLWkaCXIBK2oAAAwjqAEAMIygBgDAMMokpE0hdNIYiE5aWr7oytgDszP7/VLtuE8GPbgdrKgBADCMoAYAwDCCGgAAwyiXAAADxPfgdNb+YUUNAIBhBDUAAIZ5Cuq2tjZNnTpVw4cP18iRI3XZZZdpz549MeccO3ZMra2tqq6u1rBhwzRnzhx1dXWlddAAABQKTx31Sy+9pNbWVk2dOlUnTpzQ9773Pf393/+9du/erVNPPVWStGDBAj333HNat26dQqGQ5s2bp9mzZ+u///u/M/IDID32PXZeGt4l//dRp/ve3veMnZjW94u36A+vZfT9C8XSO66OeR6M9Ps0Ev/QWfvHU1A///zzMc8fffRRjRw5Uh0dHfrbv/1bhcNhPfLII1q7dq0+//nPS5JWrVqlCRMmaPPmzTr//PPTN3IAAApASh11OByWJI0YMUKS1NHRod7eXjU3N0fPGT9+vOrr67Vp06ZB36Onp0eRSCTmAQAAPjTkoO7v79f8+fM1ffp0nXPOOZKkzs5OlZWVqaqqKubcmpoadXZ2Dvo+bW1tCoVC0cfo0aOHOiQAAPLOkPdRt7a2ateuXdq4cWNKA1i8eLEWLlwYfR6JRAjrDEhPB43yQG/M82Vn2Z7XoXTg+d5r3/vJCUnPeffKpiyMJLcdu3RazPPyZ7b6NJL8N6Sgnjdvnp599lm9/PLLOuOMM6LHa2trdfz4cXV3d8esqru6ulRbWzvoewWDQQWDwaEMAwCAvOfpV9/OOc2bN09PPfWUXnzxRTU0NMS8PnnyZJWWlqq9vT16bM+ePTp48KCamvgbKgAAXnlaUbe2tmrt2rVav369hg8fHu2dQ6GQKioqFAqFdM0112jhwoUaMWKEKisr9Z3vfEdNTU1c8Q0AwBAEnHPupE8OBAY9vmrVKl155ZWSPrzhyb/8y7/o8ccfV09Pj2bMmKEHHnjgY3/1HS8SiSgUCukizVJJIL17VguJhU76Py/8D7+HkHHWO+pMWPD73/k9hITiO+i3v5X4t3nVKwfuSEnWURfiPup4Rb2Jo4POOrETrlcbtF7hcFiVlZUJz/W0oj6ZTC8vL9f999+v+++/38tbAwCAQXCvbwAADCOoAQAwjKAGAMAwTxeTZQMXkw2NhYvHUpUPF58V4sVl8TJ9sdmtd1+d/KQMq3i78C4mS3bxWDJcXBbLy8VkrKgBADCMoAYAwDCCGgAAw4b8oRzwVz500vG+sfGfEr6e6Q67ZdM1nr9mddMjGRhJbvN6wxEUhvgP8YhHh/3xWFEDAGAYQQ0AgGEENQAAhtFR54h87KS98tohZ6M/HjCmNbFPx87d4en9/u/27Pa5Y5YM/ECKdBvsQy8+Khc77A+qY9c4+bCvOtV90sgcVtQAABhGUAMAYBhBDQCAYdzr2yg6aamo2NR/mhkR32Fnu6M+GdnosdPJj87bWkcd3zd/UF0c87zi7b5sDmdI8n1fNff6BgAgTxDUAAAYRlADAGAY+6hx0pJ1xv19gYy+fy4o3l+e8PW+s47FPLfYSceLH6P1zjof922nKhc6aXw8VtQAABhGUAMAYBhBDQCAYXTURuTDvul86JgzLVmHjcyL77ALsbNGbmFFDQCAYQQ1AACGEdQAABhGR+2TXOik6Zy9i98nnY+ddLJ91f+3PLXOd8wt2d2nPdi+a3prWMKKGgAAwwhqAAAMI6gBADCMjhqS6KMxdLlwv/JM+6A6ds2T7c+njv/8aeQXVtQAABhGUAMAYBhBDQCAYQHnnKlyIxKJKBQK6SLNUkmg1O/hpE22903TOduQj/uok0rzX/+zva96MKnuq053Z91TGTfJGfjjXvGOrc+wLn9mq99DSKsTrlcbtF7hcFiVlZUJz2VFDQCAYQQ1AACGEdQAABjGPmoghwy4rzZ7mLNisPuBf9S7V8b+e+iPu7xmQKecqixcgvLBiOKY53531scunRbzPN8660RYUQMAYBhBDQCAYQQ1AACG0VFnSKr7pt2x4uQnJRDfJhWfeiKl98PQFMLnU2da/OdbW9hXHe+0R7112J4ZuC2C3511oM/AJPiEFTUAAIYR1AAAGEZQAwBgGB11lqTaOaeq7317/6oLoTf3u5OO33c9FJ73asff1jrF5YDFTtqr+A773Zbc3/+erLMu5E453VhRAwBgGEENAIBhBDUAAIbZKy7zhN+ddC6I783zobPOdCftx729k/XcSceU7KOYfV4uhK843/PXhH66OaXvedpqe717qr05nXTmsKIGAMAwghoAAMMIagAADKOjhhkW93p71VfjrWcv68r9nznVz8jOxX3S8b12qp21BV57c7/3gvd8aWrM8+Avt/k0ksxjRQ0AgGEENQAAhhHUAAAYlvsFGfAXe2c+lPXv+annrkvp64/HddqpdtbpuLd3qiyMIZGx22L3uv9h6sB+2eveaq/n52KnHd8Jn/Kn3L/vQa5gRQ0AgGEENQAAhhHUAAAYFnDOmbpBayQSUSgU0kWapZJAqd/DSZu9D09NfhIS8qODzrT+uBthj3/u+oTnf+ra/N0rmi7xHfRQvPLDc1MfiAeZ7qwj/+j9fubxgt19yU/KoJ4qb5+fULnW9nUAJ1yvNmi9wuGwKisrE57LihoAAMMIagAADCOoAQAwjKAGAMAwbniSJaleBFSIF6Pl48Vjybw+84GY51/+ROH9e7cg/uIurzc0Sbd0XAyWKq8XcyF9WFEDAGBYSkF91113KRAIaP78+dFjx44dU2trq6qrqzVs2DDNmTNHXV1dqY4TAICCNOSg3rZtm3784x/rM5/5TMzxBQsW6JlnntG6dev00ksv6a233tLs2bNTHigAAIVoSB31e++9p7lz52rlypW64447osfD4bAeeeQRrV27Vp///OclSatWrdKECRO0efNmnX++/z1LrorvuAuxs85F8Tc0SYZO2qai3sSvD/9ZajfXsNBBw64hrahbW1s1c+ZMNTc3xxzv6OhQb29vzPHx48ervr5emzYN/ok6PT09ikQiMQ8AAPAhzyvqJ554Qq+88oq2bRt4FXNnZ6fKyspUVVUVc7ympkadnZ2Dvl9bW5tuu+02r8MAAKAgeFpRHzp0SDfeeKPWrFmj8vLU76crSYsXL1Y4HI4+Dh06lJb3BQAgH3haUXd0dOjw4cM677zzosf6+vr08ssv69///d/1wgsv6Pjx4+ru7o5ZVXd1dam2tnbQ9wwGgwoGg0MbfQHLhc46fh/0p567LuHruchrBx2PTtp/O//t3IEHv+btPY58zVvH7NiSDA88BfXFF1+s1157LebYVVddpfHjx+umm27S6NGjVVpaqvb2ds2ZM0eStGfPHh08eFBNTU3pGzUAAAXCU1APHz5c55xzTsyxU089VdXV1dHj11xzjRYuXKgRI0aosrJS3/nOd9TU1MQV3wAADEHabyF67733qqioSHPmzFFPT49mzJihBx54IPkXAgCAAQLOOef3ID4qEokoFArpIs1SSaDU7+HkDD866gOXrMz698w3M+rO9XsIec9rf5wOdNDZV7k2tb3s2XbC9WqD1iscDquysjLhudzrGwAAwwhqAAAMI6gBADCMz6POE6l+3vXJeOGtnRn/HkCq6KSRb1hRAwBgGEENAIBhBDUAAIbRUeeIZPuks9FRA7kg1c+GHgo+Tzr7cm3fdCpYUQMAYBhBDQCAYQQ1AACG0VEbkeq9ur1+PZ02gFwWf11APnfWrKgBADCMoAYAwDCCGgAAw+iofeLH50cn+v6Dddbc2xu5yI97fQf6vJ3PvcHhBStqAAAMI6gBADCMoAYAwDA66izxu5NOZrDx9boOH0aSW/rV7+n8L3/C9n8HuWjstvK4IzsTnr/z387N1FA+Fp00UsGKGgAAwwhqAAAMI6gBADAs4Jxzfg/ioyKRiEKhkC7SLJUESv0ezpBZ76SHYu/Mh/wegnl01tk3sKNOXbp7bDrq9Mv1e3ufcL3aoPUKh8OqrKxMeC4ragAADCOoAQAwjKAGAMAw9lEDaVTk8e++z74Zu1f9kk9MTudwMET9KV4eU9SbnnEAEitqAABMI6gBADCMoAYAwDA6asAQOuv8kKzjDnjbbl+Qgt0eP+Q7j7GiBgDAMIIaAADDCGoAAAyjo8ZJ+9Rz1yV8nXuBww9/mHos5nkm7v2dbl7vUx35x/Mz+v6Z0PMl7mOfLqyoAQAwjKAGAMAwghoAAMPoqJE2yTrseIXQaU/Y8E8xz4v3e+xPb0/88pglmzyOKP/Fd9ZDckXqb5FI+IrYzjn008SdsoXOORk66cxhRQ0AgGEENQAAhhHUAAAYRkedIYHyxPepdceKszQSu7x22vGy0XHHd8xeee6kkwi42OdvLGuKeX7m0ux31vFj8MrF/VEYc0vu9+7JOudc9G5L7L/nU/50wqeRFB5W1AAAGEZQAwBgGEENAIBhBDUAAIZxMZlPkl1slgwXo6V+MVq84lNz/+KYbFxclu6Lx+L93/LY9/fj4jLrF4PFX9iVDVw85h9W1AAAGEZQAwBgGEENAIBhdNQ5Kr7jTtZZJ+tf+94vvP8U8qGTTiYdnXWmO2kM9O6VSebcJX55KNLdQbviQFrfr/yZrWl9v1zCihoAAMMIagAADCOoAQAwrPCKyTwV31kXFXsrsYbS1xZir+23+I7Za398Mp01nXTmHZ3dmPD1YKQ/4etFvd7+fAf6MlBqx0l3Jx3v2KXTYp4XUmfNihoAAMMIagAADCOoAQAwjJIxQ8Z985WY5/seOy+r37+/L3Ff5LXDHozXXptO255U+2gMLlkHnQ8y3Unjr1hRAwBgGEENAIBhBDUAAIZRGmaJ3511vPgOOx2ddTLxnbbXzroQ7s0db8yS1PZN5yI/Pn86mbe/5W3eK95OvA86F1jroAtp33Q8VtQAABhGUAMAYBhBDQCAYXTUkDT4vutM99aF2Dl7RSftD6+dtHXW+mZ4w4oaAADDCGoAAAzzHNRvvvmmrrjiClVXV6uiokITJ07U9u3bo68757R06VKNGjVKFRUVam5u1r59+9I6aAAACoWnjvrdd9/V9OnT9bnPfU6/+tWv9Dd/8zfat2+fTjvttOg599xzj+677z6tXr1aDQ0NWrJkiWbMmKHdu3ervLw87T9ArrK2rxo2BDK/nT3jLHbOH5WJ/jnb+6b7S+Pug+Dx86lzQSHvm47nKajvvvtujR49WqtWrYoea2hoiP6zc04rVqzQrbfeqlmzZkmSHnvsMdXU1Ojpp5/W5ZdfPuA9e3p61NPTE30eiUQ8/xAAAOQrT7/6/sUvfqEpU6boq1/9qkaOHKlJkyZp5cqV0dcPHDigzs5ONTc3R4+FQiE1NjZq06bB/5bd1tamUCgUfYwePXqIPwoAAPnHU1Dv379fDz74oMaNG6cXXnhB3/72t3XDDTdo9erVkqTOzk5JUk1NTczX1dTURF+Lt3jxYoXD4ejj0KFDQ/k5AADIS55+9d3f368pU6bozjvvlCRNmjRJu3bt0kMPPaSWlpYhDSAYDCoYDA7pa/OJxc7aj/uB57v4e3fHy8V902cuLbxOGsgmTyvqUaNG6dOf/nTMsQkTJujgwYOSpNraWklSV1dXzDldXV3R1wAAwMnzFNTTp0/Xnj17Yo7t3btXZ555pqQPLyyrra1Ve3t79PVIJKItW7aoqYm/1QIA4JWnX30vWLBAF1xwge688079wz/8g7Zu3aqHH35YDz/8sCQpEAho/vz5uuOOOzRu3Ljo9qy6ujpddtllmRg/AAB5zVNQT506VU899ZQWL16sZcuWqaGhQStWrNDcuXOj5yxatEjvv/++rr32WnV3d+vCCy/U888/zx5qj+isc9PYuTv8HkLGWeuk6aALY191IfP8oRyXXHKJLrnkko99PRAIaNmyZVq2bFlKAwMAANzrGwAA0whqAAAM4/OocwSddWGI738t7KtOtZOmQ86+ZJ31B9XFnt+z4u2+lMaUDPf2/nisqAEAMIygBgDAMIIaAADD6KhzFJ21DZneN52sH07WYfux53nB738XdyT2+a13X529wWRItj9/OlXxnfVQeO21M91pFxJW1AAAGEZQAwBgGEENAIBhAeecqSIxEokoFArpIs1SSaDU7+HkLAuddTK50GGvbnok4evLzrI/z9k2sKNOTaY77eqVA3v8o7MbM/o9LeqpTO+6LVlHXej7pk+4Xm3QeoXDYVVWViY8lxU1AACGEdQAABhGUAMAYBj7qPOUxX3W8eL3XXuVasedrH8eDJ109t1x008Svp4P+7ItCEZi94Z77azppDOHFTUAAIYR1AAAGEZQAwBgGB11gciFztqr+I77Py/8D59Gkt/SvS863eiw/XHao9m/j3yhYkUNAIBhBDUAAIYR1AAAGEZHXaDyobOmk5YW/eE1T+ffM3Zi0nOsd9JeJeuw7105IUsjsS1+3zQdtB2sqAEAMIygBgDAMIIaAADDCGoAAAwLOOdS+2SDNItEIgqFQrpIs1QSKPV7OAXN2gVmFi8ey/SHdHi9WCwbeh3XoCazfNGVfg8hoVOe3OL3EAreCderDVqvcDisysrKhOeyogYAwDCCGgAAwwhqAAAMo2zCx/L7pigWO+l48R1yshuKWOycvSoNnIh5TmcNZBYragAADCOoAQAwjKAGAMAwyiWctEx31rnQSR9ziff250MH7RWd9UBe9ykfnd3o6/eHbayoAQAwjKAGAMAwghoAAMMokzBk8Z11vGQddi500vCuEDvrez85IaWvp1NGIqyoAQAwjKAGAMAwghoAAMPyvzyCb5J12NqfnXGkItm+aSSXi511qp0zkE6sqAEAMIygBgDAMIIaAADD7JdFyFvLzordZ710f5JOO8Poo7PDWmdNHw3rWFEDAGAYQQ0AgGEENQAAhtFRw4xsd9Z00jZ47azplFFoWFEDAGAYQQ0AgGEENQAAhtFRAzCFDhqIxYoaAADDCGoAAAwjqAEAMIyOGmbF76uOF7/Pmn3RuemesRP9HgJgGitqAAAMI6gBADCMoAYAwDA6auSs+A570R9e82kk8IJOGvCGFTUAAIYR1AAAGEZQAwBgGB018kay7pMOOzvooIH0YkUNAIBhBDUAAIZ5Cuq+vj4tWbJEDQ0Nqqio0NixY3X77bfLORc9xzmnpUuXatSoUaqoqFBzc7P27duX9oEDAFAIPHXUd999tx588EGtXr1aZ599trZv366rrrpKoVBIN9xwgyTpnnvu0X333afVq1eroaFBS5Ys0YwZM7R7926Vl5dn5IcATkY+dtj0wUD+8xTUv/3tbzVr1izNnDlTkjRmzBg9/vjj2rp1q6QPV9MrVqzQrbfeqlmzZkmSHnvsMdXU1Ojpp5/W5ZdfnubhAwCQ3zz96vuCCy5Qe3u79u7dK0l69dVXtXHjRn3xi1+UJB04cECdnZ1qbm6Ofk0oFFJjY6M2bdo06Hv29PQoEonEPAAAwIc8rahvvvlmRSIRjR8/XsXFxerr69Py5cs1d+5cSVJnZ6ckqaamJubrampqoq/Fa2tr02233TaUsQMAkPc8BfXPf/5zrVmzRmvXrtXZZ5+tnTt3av78+aqrq1NLS8uQBrB48WItXLgw+jwSiWj06NFDei8gFfS9ACzyFNTf/e53dfPNN0e75okTJ+qNN95QW1ubWlpaVFtbK0nq6urSqFGjol/X1dWlc889d9D3DAaDCgaDQxw+AAD5zVNHffToURUVxX5JcXGx+vv7JUkNDQ2qra1Ve3t79PVIJKItW7aoqakpDcMFAKCweFpRX3rppVq+fLnq6+t19tlna8eOHfrBD36gq6++WpIUCAQ0f/583XHHHRo3blx0e1ZdXZ0uu+yyTIwfAIC85imof/SjH2nJkiW6/vrrdfjwYdXV1emf//mftXTp0ug5ixYt0vvvv69rr71W3d3duvDCC/X888+zhxoAgCEIuI/eVsyASCSiUCikizRLJYFSv4cDAEDanXC92qD1CofDqqysTHgu9/oGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMMIagAADCOoAQAwjKAGAMAwghoAAMNK/B5APOecJOmEeiXn82AAAMiAE+qV9NfMS8RcUB85ckSStFG/9HkkAABk1pEjRxQKhRKeE3AnE+dZ1N/fr7feekvOOdXX1+vQoUOqrKz0e1g5KxKJaPTo0cxjCpjD1DGHqWMO08PKPDrndOTIEdXV1amoKHELbW5FXVRUpDPOOEORSESSVFlZyX+UacA8po45TB1zmDrmMD0szGOylfT/x8VkAAAYRlADAGCY2aAOBoP613/9VwWDQb+HktOYx9Qxh6ljDlPHHKZHLs6juYvJAADAX5ldUQMAAIIaAADTCGoAAAwjqAEAMIygBgDAMLNBff/992vMmDEqLy9XY2Ojtm7d6veQzGpra9PUqVM1fPhwjRw5Updddpn27NkTc86xY8fU2tqq6upqDRs2THPmzFFXV5dPI7bvrrvuUiAQ0Pz586PHmMPk3nzzTV1xxRWqrq5WRUWFJk6cqO3bt0dfd85p6dKlGjVqlCoqKtTc3Kx9+/b5OGJb+vr6tGTJEjU0NKiiokJjx47V7bffHvPBDczhQC+//LIuvfRS1dXVKRAI6Omnn455/WTm7J133tHcuXNVWVmpqqoqXXPNNXrvvfey+FMk4Ax64oknXFlZmfvJT37i/vd//9d961vfclVVVa6rq8vvoZk0Y8YMt2rVKrdr1y63c+dO96UvfcnV19e79957L3rOdddd50aPHu3a29vd9u3b3fnnn+8uuOACH0dt19atW92YMWPcZz7zGXfjjTdGjzOHib3zzjvuzDPPdFdeeaXbsmWL279/v3vhhRfc73//++g5d911lwuFQu7pp592r776qvvyl7/sGhoa3AcffODjyO1Yvny5q66uds8++6w7cOCAW7dunRs2bJj74Q9/GD2HORzol7/8pbvlllvck08+6SS5p556Kub1k5mzL3zhC+6zn/2s27x5s/uv//ov98lPftJ9/etfz/JPMjiTQT1t2jTX2toafd7X1+fq6upcW1ubj6PKHYcPH3aS3EsvveScc667u9uVlpa6devWRc/53e9+5yS5TZs2+TVMk44cOeLGjRvnfv3rX7u/+7u/iwY1c5jcTTfd5C688MKPfb2/v9/V1ta673//+9Fj3d3dLhgMuscffzwbQzRv5syZ7uqrr445Nnv2bDd37lznHHN4MuKD+mTmbPfu3U6S27ZtW/ScX/3qVy4QCLg333wza2P/OOZ+9X38+HF1dHSoubk5eqyoqEjNzc3atGmTjyPLHeFwWJI0YsQISVJHR4d6e3tj5nT8+PGqr69nTuO0trZq5syZMXMlMYcn4xe/+IWmTJmir371qxo5cqQmTZqklStXRl8/cOCAOjs7Y+YwFAqpsbGROfyLCy64QO3t7dq7d68k6dVXX9XGjRv1xS9+URJzOBQnM2ebNm1SVVWVpkyZEj2nublZRUVF2rJlS9bHHM/cp2f9+c9/Vl9fn2pqamKO19TU6PXXX/dpVLmjv79f8+fP1/Tp03XOOedIkjo7O1VWVqaqqqqYc2tqatTZ2enDKG164okn9Morr2jbtm0DXmMOk9u/f78efPBBLVy4UN/73ve0bds23XDDDSorK1NLS0t0ngb7s80cfujmm29WJBLR+PHjVVxcrL6+Pi1fvlxz586VJOZwCE5mzjo7OzVy5MiY10tKSjRixAgT82ouqJGa1tZW7dq1Sxs3bvR7KDnl0KFDuvHGG/XrX/9a5eXlfg8nJ/X392vKlCm68847JUmTJk3Srl279NBDD6mlpcXn0eWGn//851qzZo3Wrl2rs88+Wzt37tT8+fNVV1fHHBYwc7/6Pv3001VcXDzgatquri7V1tb6NKrcMG/ePD377LP6zW9+ozPOOCN6vLa2VsePH1d3d3fM+czpX3V0dOjw4cM677zzVFJSopKSEr300ku67777VFJSopqaGuYwiVGjRunTn/50zLEJEybo4MGDkhSdJ/5sf7zvfve7uvnmm3X55Zdr4sSJ+sY3vqEFCxaora1NEnM4FCczZ7W1tTp8+HDM6ydOnNA777xjYl7NBXVZWZkmT56s9vb26LH+/n61t7erqanJx5HZ5ZzTvHnz9NRTT+nFF19UQ0NDzOuTJ09WaWlpzJzu2bNHBw8eZE7/4uKLL9Zrr72mnTt3Rh9TpkzR3Llzo//MHCY2ffr0AdsC9+7dqzPPPFOS1NDQoNra2pg5jEQi2rJlC3P4F0ePHlVRUez/louLi9Xf3y+JORyKk5mzpqYmdXd3q6OjI3rOiy++qP7+fjU2NmZ9zAP4fTXbYJ544gkXDAbdo48+6nbv3u2uvfZaV1VV5To7O/0emknf/va3XSgUchs2bHB//OMfo4+jR49Gz7nuuutcfX29e/HFF9327dtdU1OTa2pq8nHU9n30qm/nmMNktm7d6kpKStzy5cvdvn373Jo1a9wpp5zifvrTn0bPueuuu1xVVZVbv369+5//+R83a9asgt9a9FEtLS3uE5/4RHR71pNPPulOP/10t2jRoug5zOFAR44ccTt27HA7duxwktwPfvADt2PHDvfGG284505uzr7whS+4SZMmuS1btriNGze6cePGsT0rmR/96Eeuvr7elZWVuWnTprnNmzf7PSSzJA36WLVqVfScDz74wF1//fXutNNOc6eccor7yle+4v74xz/6N+gcEB/UzGFyzzzzjDvnnHNcMBh048ePdw8//HDM6/39/W7JkiWupqbGBYNBd/HFF7s9e/b4NFp7IpGIu/HGG119fb0rLy93Z511lrvllltcT09P9BzmcKDf/OY3g/4/sKWlxTl3cnP29ttvu69//etu2LBhrrKy0l111VXuyJEjPvw0A/F51AAAGGauowYAAH9FUAMAYBhBDQCAYQQ1AACGEdQAABhGUAMAYBhBDQCAYQQ1AACGEdQAABhGUAMAYBhBDQCAYf8P2s7CGmKzovoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "atlasarr = atlas.get_fdata()[:,:,30]\n",
    "plt.imshow(atlasarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The content of such array can be summarized to the following numbers, each of which is associated to a brain's zone, as shown in the [original web page for this atlas](https://neurovault.org/images/14255/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0. 2101. 2102. 2201. 2202. 2321. 2322. 2501. 2502. 2611. 2612. 2701.\n",
      " 2702. 2811. 2812. 2821. 2822. 2831. 2832. 3001. 3002. 4101. 4102. 4111.\n",
      " 4112. 4201. 4202. 5001. 5021. 5022. 5301. 5302. 5401. 5402. 7001. 7002.\n",
      " 8111. 8112. 8121. 8122. 8201. 8202. 8301. 8302. 9021. 9022. 9031. 9032.\n",
      " 9041. 9042. 9110. 9120. 9130.]\n"
     ]
    }
   ],
   "source": [
    "print(np.unique(atlasarr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we can simply mask a zone by searching for it in the array. Here, I'm using numpy's *masked arrays* to do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x12cba96f0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAADsCAYAAADUxMRPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOyUlEQVR4nO39eZxcVZ34/7/OvbV3Ve/pLenudBayQcISCAEFlAwZZGZAccfPIPoZXCIKcUbl9xEQB2VEx3FUlMFxUH8qLuMAiqOIYRPIAklYQshKyN7d6bV6q+Xee75/VKfSnV6ru7bufj991EPq1q17z63qnPuu99mU1lojhBBCCJElRq4LIIQQQoiZRYIPIYQQQmSVBB9CCCGEyCoJPoQQQgiRVRJ8CCGEECKrJPgQQgghRFZJ8CGEEEKIrJLgQwghhBBZJcGHEEIIIbJKgg8hhBBCZFXGgo97772XuXPn4vP5WLVqFVu2bMnUqYQQ04TUG0LMDBkJPn75y1+yfv167rjjDrZt28aKFStYu3Ytzc3NmTidEGIakHpDiJlDZWJhuVWrVnH++efz3e9+FwDHcaitreWmm27iC1/4wqjvdRyHY8eOEQqFUEqlu2hCiHHQWtPV1UVNTQ2GkZ3W2cnUGyf3l7pDiNxJpd5wpfvksViMrVu3cuuttya3GYbBmjVr2Lhx45D9o9Eo0Wg0+fzo0aMsXbo03cUSQkzA4cOHmTNnTsbPk2q9AVJ3CJGvxlNvpD34aGlpwbZtKisrB22vrKxk165dQ/a/++67ufPOO4dsfwvvwIU73cUTQoyDRZxn+V9CoVBWzpdqvQFSdwiRb1KpN9IefKTq1ltvZf369cnn4XCY2tpaXLhxKalAhMiJ/sbYfG6+kLpDiDyTQr2R9uCjvLwc0zRpamoatL2pqYmqqqoh+3u9Xrxeb7qLIYSYQlKtN0DqDiGmsrT3JPN4PJx33nls2LAhuc1xHDZs2MDq1avTfTohxDQg9YYQM0tGml3Wr1/P9ddfz8qVK7ngggv41re+RU9PDzfccEMmTieEmAak3hBi5shI8PG+972PEydOcPvtt9PY2MjZZ5/NH//4xyGdyYQQ4iSpN4SYOTIyz8dkhMNhioqKuIyrpdOYEDli6ThP8QidnZ0UFhbmujjjInWHELmVSr0ha7sIIYQQIqsk+BBCCCFEVknwIYQQQoiskuBDCCGEEFklwYcQQgghskqCDyGEEEJklQQfQgghhMgqCT6EEEIIkVUSfAghhBAiqzIyvbrIb0YgQM8VZ9JbbiY25GCOW+XArE0t2K/vzf7JhRBC5JQEHzOQUVJM/GOt3Lnw98S1C1tnPwHWYQf4/jfeSZkEH0IIMeNI8DENuaoqCV80F9ujhn09Wqg4u+RVQkYERxvYDL9fJrmVRftZGs/7LhzyWuGeLvT217JeJiGEENkhwcc01HfmHC66bTPLA4eHfd3EodjsxclBxuOkAiPK3e94kN6/9hLTLhx9KgD61n//HfXbc1Y0IYQQGSbBxzTiqqqkd0UtLSs8LPEfY677xIj7xrUrJxmPgcrMbsrM7mSzj9Pf/zlWHyX6jvOTfVH8x3twXn4d8msBZiGEEBMkwcc00nVhPavv2MIy/xHmelpyXZxxM5WT+A+dCEC+tOq3HD+3JPn6fzzzdhZ91osTieSohEIIIdJJgo8pzDW7hti8yuSA6bbFLs4rOMBsV3tuCzZJFa4uKlxdyecFNV3ELl6GEXMG7Wd2x2DnPnQ0mu0iCiGEmAQJPqawpivruXzdRgJGDIBydxdVrs4cl2riTOVg4mBrI9kEA/DZJX/m9W/UDOoXAvD7N5bR8NkKrIPD920RQgiRnyT4mEJcs2uwK041R3TVwxWFr+JT8RyWKvOqXJ1UhYYGVQerSmleOg9/SWjQdhW30W8ewenpyVYRhRBCpECCjynk0HVzufQ9W5PPr/CfmPaBx2iunvUSG7/cTZ/tHrR9f7gcz5cWoJ5/OUclE0IIMRoJPvKYq6oSXRhMPu+ps/lA2cYclii/lJnd/E3JS0O2bwvM5bcNb6P0xLyhb2pqwQ6HM184IYQQI5LgI08pt4d9n5rHWZeemgH0E8PcaKeT0/t6TNQi33GW3bSDpsjg5pi4bdJ1/1JCv9w06XMIIYSYOAk+8oVSmGWlKK838dzjxp7Xxy2z/5Tbck1BBUaU95S/MGR7j+PlS7UfomR2TXKbE+7C6eoasq8QQojMkeAjT5hFhez+/y2k7qzjABhK88mqZ3Jcquw4fZKxTPEZMS549yvsWlOR3Bb5n3rK75emLCGEyCYJPnJNKcxQCCrKqVrazF3zH8p1iXIi04EHgIlOZETKT227peajGT+vEEKIwST4yDGzYhZ7PjePkkVt/H3t5lwXRwghhMg4CT5ySHm9UFzIGece4v/VP5rr4oiTDHPk17Qja8wIIcQkSfCRI645s9nz6Trc87pYV5m9jIeNkdPVbAfKVl+PVJiVFdj1lWAMv+ie0RtD7z4gU7oLIcQkSPCRC4aJUxpi1Vte5xNVT2T11I42cr6a7UC5CjxsVCL40YBSp7IZoQL6avzoEYIPT6cbz36XBB9CCDEJEnxkmWtuHXtvnI1dF+GfiveO/YZpKB8yHnfvewfH9syi9JhGnbs0uT0W8qLVyMGZ7TNwls1D2c6I+wzigHGkGfvEickWWQghpg0JPrLMnlXEO9a+wHtKtuS6KDmVy8DDRnH0QDlznkqUpKc+OMY7BrzXa9A7JzD+k2kIdfWCBB9CCJGU0h3g7rvv5vzzzycUClFRUcE111zD7t27B+0TiURYt24dZWVlBINBrr32WpqamtJaaCGmknh1McbyxZiVFWPvPE1J3SGEGCil4OPpp59m3bp1bNq0iccff5x4PM4VV1xBz4DVQ2+55RZ+97vf8etf/5qnn36aY8eO8a53vSvtBRdiSlAQmeWhe34RlBblujQ5I3WHEGIgpfXExw2eOHGCiooKnn76aS655BI6OzuZNWsWP//5z3n3u98NwK5du1iyZAkbN27kwgsvHPOY4XCYoqIiLuNqXMo95v75xjVvLh3nVY7YYbF7jsF1H36ctxTsHvb1TLH748zROpza2uAfd7ybyK7i8R93ToS7L/gfyszu8b8nTWu4jFdMm3zv8Ns4cKIM2zbQtsK/w0/JXjtrZUCD70QUd2sPtHViNzVn79wTYOk4T/EInZ2dFBYWpv34UncIMf2kUm9Mqs9HZ2cnAKWlpQBs3bqVeDzOmjVrkvssXryYurq6ESuQaDRKdMDIgfAUX3G0b0E5x660Ua7hOyQWhCLM8bRluVTjG+USx8T8cwkN331+3Mftu+YC2s4NphR8ZFtEu9n7ci1lLytcfRozroEsBh6QyIBUeInM8hLca0KeBx+ZJnWHEDPbhIMPx3G4+eabufjiiznzzDMBaGxsxOPxUFxcPGjfyspKGhsbhz3O3XffzZ133jnRYuScuWgB3YtLk88757pQxsjDMJXKzQRVhnI4FCvny9v+BqfNO/xODjS8ntoQ0oI3u/nG7/8Ox39asBWKc8vKDSz1HZlgiVP3Zqycr2+9At3uSWwwwVXWh8dj420x8mNyMAVWsQ/3ggboCGO3tOa6RFkndYcQYsLBx7p169ixYwfPPvvspApw6623sn79+uTzcDhMbW3tpI6ZTZ0ryjn+jvipDcoij6bRSDJx2B+tpOE/FMZfXhx5Rz3OIaT9nJd2Mv+VoTOCGmedwebvN2Q1+NjZW8P87zuoTYnrM3xe+i5bRk+lC5UPgUe/aLmHWPEsAoe8MAODD6k7hBATCj4+9alP8eijj/LMM88wZ86c5PaqqipisRgdHR2DfsE0NTVRVVU17LG8Xi9e7wi/xPOYecZ8InXFdNUaYIz/xhaLuXiqYzFzfJWsKthPyOjLYCmHUpYGJ81NDsMcz2gNs/mZJTxfPQ8A0+1w/ZmbuCCwP239PV6NzOEHr74F2zaoKAvT3hWgqshDYH49+mgjOm6hHMhRsml0CuygF3d9Lbq7B7s1+01xuSB1hxACUgw+tNbcdNNNPPTQQzz11FM0NDQMev28887D7XazYcMGrr32WgB2797NoUOHWL16dfpKnQdaV1XQ/Lb4qE0sw4lHXTyx9wz8gRhzlrSx2Jvd4CNbrCNHmX9HC6p/wi6jcha//d5ZrFx0IG3n2NC8mAVf7kV19dLytjoCBdA9GyKlFZQ9FcFubknbuTIhWuYlWlqN/2gPtLXnR7NQhkjdIYQYKKXgY926dfz85z/nkUceIRQKJdtii4qK8Pv9FBUV8dGPfpT169dTWlpKYWEhN910E6tXrx5Xb/WpwFzQQLyqiN4qNWKn0tForUBDPG6ytXsuR2Klg143lcOZ/iPMMtPfea7c1c2xtwYoqVxF4fNvYjVmdg4FHY1y8naqOsN0bpvDzT3vBcAwNFfNe43LC3eO+3h7o1U82XoGlk408+w9VsH8YgfTk2hWMWwFChyXwq4uxQgGsAIG+mQzmCK3TWL6tP+HRHkDbtyza9A9vdjt7bkoWcZJ3SGEGCilobZqhGmnH3jgAT784Q8DiYmCPvvZz/Lggw8SjUZZu3Yt3/ve90ZMnZ4ur4fLKUXrRy6k9UILXA7KnPgvVaX0sO83TYcPLn6RNaEdkynpsGwMTliFbOmZxwv/dB7uP29N+zlGYxQUoNyJeFf5fLx+z2y+d9HPxv3++49dwr5HFmKeTDY54Ipo0KBNTgUWGlR/S5A2GBRwuHs0rmhuMgzK0clyDaITr/ma+9DbXk9/s9gEpHuo7YyvO4SYATI21HY8cYrP5+Pee+/l3nvvTeXQ+U0pXHPrsMtCRMoVyjv5m4PWCm0NUyFrxd6eCkLmguSmek8LVa6OcR+7y/GzN1o1ZGitW9ks9DQyx9POFlf2pzd3Bkwopdx9+HbP50vlf0t5oIdSb88o70x4vakKXx/9Q2UT9NC+rqBAD/eXna+tGgq0qbB9LjwV5eiT68ZYFnZHx7RojpmxdYcQYliytss4KJeb42traD/bBncso5l7x1ZsfnMuW4y65LZ3LNyZ0low+2MV/Ofei4jGBn+9AV+MdQufTltZJ0PHY8y99zXUf/kJX1DH4Yax/xR91uDAY7qJFXuIn1OX7CDrDscwtvbKCrpCiGlHgo/RKIVrzmyckiDRUoUKWFk5rW0ZnJz5XinNkd5iXvUPHUJY6upmtmtoHwFHG8TiLqz44K83Ymj2RyqIOi56Z7nwL1qAPtaE09U15BhmSQmqKJR64eMWdlMz2hr7s7I7OqGjk8CRMhzXqcXdLJ/CCmQmxNNGok+IsnXmR8GcfvyxzqdAu9SpriFuA6VU3iZshBBioiT4GIXh9XLs6jo6F9lof2YzHiPRWvHK0Rpeb64c8tqKqmN8vPpJTMbX8TUWc/HQvhUARN5q0Xp2OQse9MGLQ/uX9F0wn5YVnpTL623TVDwUSW3o6I59hPacOpd91jzalgbS3zlUJQIb2wuuPjLe90NpPSjgUKn3TxZCiGlJgo/hKIWrsgJdUkikFCiK52yQhFIa2zJx7KGdG05EgrwRq8A4Lfg4HisavpuAVsT7m2LMAgvbpemrDhBsqO9/XeO0tuN0dWH7DGKh1G/OylJQWY7LkwgmtOPgtLaNmgnR0eigpgUzHMXT5R/zXNqEuF+ltjyi0R8PpPML7f+YTnZuVU7/Nj25gEMbChUKYZgmTm/vtOj7IYQQIMHHsIxgkKPvnU/XXAcnlMPAw9CMMEgAgEMtJXyz6/Ih223bwLKG64k5mOGxOXSlgfH26sT5HJj7aDnmk9smXOZ4CA6/oxzllAPg6tVU/96DdXj8M53qvQcoPhoYe8fyElovrMT2TbS0adA/UgUFltfAMRXuXgfD1pPOdFgBE2dFHa4eC/PlvYkARAghpgEJPoahTINoMejy9DW1nFzTRevxHXE8a8DYlknvOIKM4Th2Yr4RVWChCxKBjgb6yj0UlZdheyZ25Y5LExuwcrztVTilIczuEpyurnH1BdHRKPY4OlmaLheenllYdqKs2gDLm2ImJF1O+7rS0cSiDYXtNzFsjTlaFCqEEFOMBB8ZppQelOJXSqOd0W8kqn+69kzdb7QDqs2DEe0/gQKrxMLwW7QsN+ipWkQ8COkYm+p44ejlJbi7i6l6/BjWgYOTPmby2O3thJ7ZB2YiAFMFAdovqCIWys2N2tWnQWmULc0jQggxGgk+BlIKIxhEFRcNP3/ExA972gY9YgbkZMYjkz90tVYYcYURO3lSwFZoW2EFNemc8N0xNdESsH0K7Uu9A+totGUNWhXWCIVw9VVge0xsN1nNgCinv/lFCCHEmCT4GMAsLub4B5bQW62JVWSwr4cavVkl6xl2Da4OE8ImKjujiTPC6ekluOlNgqECOs+tJFokTRVCCJGPJPiARMbD60UVheip1VizoxntZJoXzfdKg6E4OVDGiGf+lNrjwvD5cGJxcGyU24MyR05PaMsaVx+RJMfGbmrG6ArgWjyLeMDEMclNHxAhhBAjkuADMMvLabx2AZFZinhpejIep/f1yCfK0NglFralcHW4TjW/ZJDjgcaLS/AsL6b8ueNYbx4mfslZdM4duSmmZF8E45mXUh5i6kSiBLcfoaDAT3h5OZESiT6EECKfSPABKJ+X8HxwKiNpjRfyIsMxDKXA9NloB3TXwBXZMscxNb3VEC1TlOwswDjmpqfaTde8kd/jDXsIjvzyKCezsY4eQ3m9mGeUTbTIQgghMkSCjzTL54xHPnBc0Hx+Ie6l59BTo8jVam+WV+G4FGZMY8YmUAZ9agZTmblUCCFSI8GHUjBKv4OJHlIMTxuantn9T8aay0QplGmibXvCs3sqrRP9WhSDgkLtAtsFyoYJD2zKZuBhSNOREGL6mNHBh1lZQesV84mUK+yQlVLCQjIcmdc1x0Bfu5Lg4T7Uph3g2Cm9X8ctCnY24z9aQNeC0LB9P2yPQhsKMz6+DIhWifegEk/GszidGdMYk1iN1/KZuBbVY/bG0G8cwolEJnwsIYTIBzM6+KAwyImVGkomMLpFzawMx8Cb7DgnaZ20aJkmWqZw3H6Ktyh0qlkGx8Y6cBDlcuGuOodI8YDgo/8akhkQDeZ4Ot4qcNwqsY7LOOMJZetJjSbSLkVfdQCzz4vvmBck+BBCTHEzO/gQ4xY4DsV7Y/RUe+hcoNDG1JlQS9s2/n0teE/0rxejFL2zA/SVmYlAygDHVFjecRxM9QdfKVz+wOzKZDIgjtvAmT8Hs3sW+vAxWetFCDFlSfAhxqWgycH9560UrV5BeH5BjrqJTpDW2PsODNrk86wgUhJIrnKrzcQ6NOM7Xmqnd1wKpz+7MpngQ7sUfTUFGDE//pY2kOBDCDFFzcjgw6ysoOuiBnpnGWivPWyTi7bVqO0LWvWvOGtoDHNK3YpPUWAHHRyPwuw1sjLR2ET0lRv4/uocfM296O27Uu77MRz3sXZKLAeMRFYiVuKht9w1dj+eSXzVthvAwLAmkQHRY/fTFUKIfDcjgw9dWcrRNUBolL4etoJRF4BTifuQ24EpGnwoBWYojtbgWF6MeH52YonM0hyf5aJwfyFlO1zo6OSDD+vNQ/Bm4r8V4F+xhN7y4tHfNMmv+WQGxBWZeAZEaVm4Tggx9c3I4AP6Mxf0ZziG32F8B3IUTjwNwyCVRpk6q51YtQan1wWWwjXBwMPTpSh80xnxhui4FZ0NBlZghBumZsxsg6ddUXjIIXA8CvbkA4/hGK1hivd5iRW56Z01TAYkjfd7xwWWT2HYA4IQiSeEEDPIjA0+gESFbxmTq/j7V4OdNAV4bVQ2sygazC4Ts3fi5fe2aQp//+qInR/N0hL63r8YKzD8+ZNGKUJBo0PooW3oeCxj92jryFHUkaMElp5BX3lZRmOBRAZE4Yr2N78MmLBMCCFmgpkZfDigbIW282viJm0rHE3GMyBag9PnAlvhssZ3or5Sg8CFy+mp9Q9NCmk94iRgui9C0RtxfO2Jqby0AT01BrHQafuPFYikMs7WMHHVVKH9pw1f0RqaW7HD4RHfqjq7CR0KTnoOF60UkTIX8cAYf2P6tP8XQogZYEYGH0prlKXQ47zxZkV/FkZDxjMg2lGYHS7MKOO+6XXXKnrmBNGKlIbZOr29eB/bhs9IfNbK68V5z3JioYEFGuaNk/hqDI+b3jNriJSeNnephpKtBowSfFhHj2Eeb5r4yfsp04S3nkU8MPLCeZLxEELMVDMy+ADys8Lv7/+gbUVGZ+12FNqtcYbp16KsxJTjQ4pm6Il/ZI49KHERPBZHKzd9FYp4cOhRPR0KX9up7QXH42gnhbMbBtpIZFkG0WCVBXHPrUO3dYycAUnDaBqtHTytvQRdQz/jWMgkFpx41k2bCkqLMV0unLYOdDwLyxILIUQazdzgI1+dzIBYGT6N38bxD91uhl2YmenTmTivZeF54iXKvV7arzmL9sVDb84l+ywK/vgK+mRTjm2j0xAQoKCrwY+q91P8ig92jpwBmTStcV7dg+e1065PGXjOXkTb0gmt1wsk+oz0zS/DiJXgfVVjnzgxycIKIUR2SfCRj8b4ke8Kxikp6kn5sHHLpLO9AB0zRmzW0G6NYyuUM3wGJFmGXoW3XRM4YScWfkuBtiy0owmcsIgXuIe87jsRm9D6JYbPhz5rIfGgG8PS+FptYoUmzoBTnMyIWMV+3LNrhpatpwe7ozPlcw/rtIyPUVCA8vuwXEZittNJDJnVRmLSMSGEmIok+JiCzqs/xD21vyXVxP1BK8DHXvo/9LQMN/QkwfHbOD6F2WOMOgom0KiZ9dvd6J5enGg0xZIAjo3n6VepeG7on6COxSfUxGNUVfD6p9yEivqo+IYX3ytvoC5aRF/50HVru+r9GHPqhmwPHuqFza9OeBXd0ag51UTqi3FMhSfsSH8PIcSMNSODD9UXJdBoECtSxMpsyOd1ShS4gzH8/lPt+gsKTlBpejFSDD/iupe6knaOGqP3KNFAV3MQI37qz8PVpfAOSAgEWiyczjDamnj7kI5G0RMJXPoZBQWwsB7tTgQX3ZV+Sko7qAx209pQT7E1F2uE0SbaBHtoTIIVcOOtmAV9kVFHxQyhFGZRIbhH7mDqBL04/dkKlUofluQ5QBunAkIHoDiEqR2czi7p+yGEmDJmZPBhHzpC7X914cytZt/7g9iFGezkkAZXLNjFpyueSD4vNsDAl/JxKk0P983/FZExJlCzteLWg+/k5dfrk9sK95lUPLwnmRHQkSjOJAKPtFhYz/7Pu6kqTQQJphFmeagNl3KwP3yYlt4A0eeCBI+Mv/tub6WbaGkDgWMR1KYd4+58qjwe4ivmESkdOfgwrMnNTqqVwvIbgzrSWstnYcbKCW4/inXk6ISPLYQQ2TQjgw9tWdgtrbgKAvhaConFRrgZK7BCDtqT0bEnI57bLIjj8VgsCDTR4Eo92DidW5lUm8P0Mj2Ng8PCUDO7yyqItPlwt7rwdTjYLa2TLkM6OR6TqtIwZ5YeH/La/MIWSnwBXrPLcfdobG9iavOxaBMsU+F4TQYlRgwTM1gAxgh/K14vkYALyz9yYGdGwZxEcuJkf5WBwYftVYnRL65h0jhCCJGnJjXL1r/8y7+glOLmm29ObotEIqxbt46ysjKCwSDXXnstTU2TnzchE+yjx6n/8Rss/P7hYR8L/quJwJEcVepKc/WiV/j5eT/kg4WvZf30N5X/hZ+u/CEqarDw+4cp3LAn62WYrIjlpmSvReiZvfjaJpelMUuL6btoEd2XjvBYPZdY4eh/K7ZHES8wJvyw/Ma4Z/3PZ1O93hBCTN6EMx8vvPAC//Ef/8Hy5csHbb/lllv4/e9/z69//WuKior41Kc+xbve9S6ee+65SRc23bRlYR1vHPF1w+fD21ZJpMyFE3DQ7ixmQBRUeTtZ5nGRiwRVm+1mf3wWnnYT6/CRrJ8/Xcw+B7u9EyM++Lsz4olmkJEoS2MGC5IjeVQoSDxoYHsmcfcffaHkGWE61BtCiMmb0F2tu7ub6667jh/84Afcddddye2dnZ388Ic/5Oc//zlvf/vbAXjggQdYsmQJmzZt4sILL0xPqbPEiUap+v1BKouCHL6qnO6FOWh+yTIHB1tr3rlhHfN/6jD/0DFy3LMjLZSlMU5eiIbC/d0Ybw5trkkqK6H34kXJIcmOWyU7i4qJmSn1hhBibBNqdlm3bh1XXXUVa9asGbR969atxOPxQdsXL15MXV0dGzduHPZY0WiUcDg86JE3tMY6egxnzxt42zRGrwnpmpLd0CifjeG3EqNtFCivjRGwMAIWpt8iYGRn9IIz4H8HrRhboj4C+z2YT23DeuPNrJQhnWytaIsFaIsEQCkMvy/xl67pn9IcVMzC6eoe8tDdPejuHtCaWKFJtDjxiBcYQ2dMFSlJZ70BeV53CCFGlXLm4xe/+AXbtm3jhRdeGPJaY2MjHo+H4uLiQdsrKytpbBy+eePuu+/mzjvvTLUYWaUti8o/HWbWiyGOvb2E8JL4pI9pBiw+tvwvlJo9fHPn5fR1+XjP8m1cEtqV3OcsTwswdgfRyXAGTOQe0RZr/3ALdb+H+t1N5PcYoJG1xwJsf/4M/I2KwgIbtXox8dCpP3WtoLeuELPizKFvNhRageORYCOd0l1vwNSoO4QQw0sp+Dh8+DCf+cxnePzxx/H5Jj/6AuDWW29l/fr1yefhcJja2tq0HDudrMNH4IjCvXKCKWBDo1z9P7sBny/OpQW7qDKj/GfgYuIxF6uD+7jCP3Dm0swGHsCguUIcrQnud+H73fNTIvBQjqazz0djXwiPadO/LB+tkQL8jYrgUQfHpYiUDZ1F9WQnzqEHHWZNmHylQTn9fVcyMClaumSi3oCpU3cIIYZKKfjYunUrzc3NnHvuuclttm3zzDPP8N3vfpfHHnuMWCxGR0fHoF8xTU1NVFVVDXtMr9eL1+sd9rXpxFcU5ZNLn2GWK5EaDpl9LHJbuPFw1xkP0+X4ucjXhFsVDHqfncpS8pPkVlPlrptg7D9K+bfn0Vwyj855Jnb/FBvKgYLWcd6Mp3AnUMOGwt2dGK1hnDwbBj1QJuoNmDl1hxDTUUrBx+WXX86rr746aNsNN9zA4sWL+fznP09tbS1ut5sNGzZw7bXXArB7924OHTrE6tWr01fqHDJswE7hbqUAt0PQH+WdodeY4xq4oFgis3G53wa6gYIhbzezGBAYeuQ1X/KR3d6O68lteEtLMKMLhs9kjMcUuuZBtMZo6cQ6eizXJRmV1BtCiNOlFHyEQiHOPHNwO3lBQQFlZWXJ7R/96EdZv349paWlFBYWctNNN7F69erp0WNda8qfb6bwzaJxv6VluY933PAsFwTfoNQYefZLkTojFEIvnkuswI3tTTGCmMIZj6lmxtcbQogh0j6BxL/9279hGAbXXnst0WiUtWvX8r3vfS/dp8kepWBA9sHesx8zhfm2guWreG/xC5zt9QISfKST8nrorfJj+RPfj0qh24NWTN2MByRH7UwX067eEEKMatLBx1NPPTXouc/n49577+Xee++d7KFzTrk9RNasoLsm8TEpDeUvtuO8smuMd4ps0N09BF9rRk9gavFIXTF9FUM7ok4FRlwTOtCD0dmL09E59hvy0HSuN4QQY5uRa7uMl3K7aFvipmtx/9BaR1HQGML7Sm7LJRKcSARngvOQuIvPmrrBhw3G4WbspuZcF0UIISZEgo8BjIIC+i5bSqQk8UvaMRWRigG5baVpOctNsHRoO3Th/j7UxpezVVQxgzkuRXxhDa7qctQbR7Blci0hxBQjwccAKhCgaaWbyOwRJhFT0DM/Ts/8oS/ZngClm1Rez7cgpgfHBd11Poy4l+ITQZDgQwgxxUjwQWLUROSiRfSVu4gHJxY89MxWeK+9gIIjfbD5VdS5SzlwTSFqUTc1LguQ+QjSzQgEULU1aHfqf8bRYi9qmClUpnxHVCGEmAIk+ACMokKOvdVNbJY14RtPX22co3Og5KUCKl40aV1eyOPXf506V5Dh5u8Qk6cCfroXl2J7U5/fQytgmDhTIUNwhRAi02Z08GEWFxE9dwFds9zYfj35X7wK+mYp+q48l85F4Ca7M5ROVJ+O8aH9V/PyG3Oo3zH5dWsyzVU7h5bLaokVKXpqNLr/r9iIQeF+cPdOsOlL5v4QQoismNHBhyot4dAVHqxiO22TJkTmxDk0G+objuNTBs5wP69zbOBicgAdjsWhX8xj0Q9eQNv5v6pLdEEFhR8+wtxg26Dtb3aXcqK5dmLBx1Ra00UIIaa4GRl8mMVFWMsa6Kz24nh1Ykn7dFLQ1uvn/o4VLPA2sTbQjFflz7BOA2NQAOJTio5VMbS6gKrnO3Be2pnD0o2P27BxG4MDpZA7wv75Dn2zEukLw1YUHNZ4uvMvAJwMZYO/JY7Za6F7+3JdHCGESNmMDD6oqeTA1X6sQpthex2mQdeJIPe1XkrN7DYuWvoTvGb+BB8wOAAJGR62Xv4d2t7m8M5//xzVL+W2bBM1y9fN21btwOlvP+uM+dj9h4XTLvgw4xrfa0ewmpqxZXSVEGIKmpmJZqVwXP0Zj0y18WvAVthO/n7ExoD/BQ0vlaaL7hUROv5+NcaKJbkuXspMpfGaFn4zjt+M4zHtxPc72mOK0lrLsG4hxJQ1MzMfYlhuZfLc275Nx6UG1/7wH6mdBnOmaTMxL8ZIVOaSX0IIIUaQvz/LRdYZGJSbfua53UQX9dH7zlWYyxbluliZ1Z8ByetRLjoRJBlW4qFsDY5kPYQQU5cEH2IIA4M/v/U7fP1fv8e+D5XmujgZpY1EZkSbif/OxyBEOeDudXCdfPTY4OT/qCQhhBiJNLtkWF/MzfORGsrMbgACRpSlbjuvRr8MZ47LzxwXiflP8oyrO87rb9RwtKyIFZVHCbpiaTu2GdO4e4a2w2hTES8whgzHVQ64e5xENmKCHLfC9hojBj7K0YnsR3+2I02jwoUQImck+MiwcFsBn3/xXcm+jWXF3Tyw9CcsyO/YI7+9soeld5TTs6KGo7cUs6gofau7+k/EcW/dC87gAMSYVUbnympsz+AIwYxpgi8dw2lrn/A5VU0lvWeUodXIaRclnUuFENOIBB+ZZins7lORRofLz6uxanr1iUG7mWjqXTZBI3trwDg4HLGidDluItpFDHO4nXBVVyWf6r4+7I7OrJVxODoaxTpyFH9piN1HZtHSPfz09XHLxIykeHAFyuNGR6I4PT2nNvt9uLsdTPdpwUfUQXf3DNp33Kdye1C+/u/bASUpDSHEDCHBR5ZFu7x8cdvVGKdNbOb1WHztzP/hcn9v1srS68S56/iVvNQ0m872AnTf0ODD06c4/IF5yecleyx8/7s1P/oc7HmTxXdVgzlM0ATgcuheYBELjb9rU2+lh+glCwkc74MtryWv027rwP/8bjAGBx/adrB7JvadmVUVROdX4LjUlB72K4QQqZq5wcfp8zxk60enrYiHh2Y3Yl6bXdFq6l0jzy5aZmpKDN+4TxXVcU7YFvER7mxdjpv9neV0tBVgtrkxe4ffL1py6sOJlJgUzCpD9/bhdHWNuywDGaEQyjNKu1N/84Pui4yaUXAiEdh3YOTDuD24Zi8fNfgwLDDig798bYA21OBPzbGxJ7l0vfJ6UR5PoklHa3TAhxUYIXASQohpbEYGH9pUaLdGuU+16+u4kb0AZLgyxQ2+v+MS/su7esR9rpv/IutLd437mG9Y8IUD76U94h/2dcs2aTlcjCtsYljjO2Z3rSL63gUUHrII/H4b2hrnG/spt4eey5cQrh1t8o3EqJPSXTE8f9qa0cm0gof6cO05MmS7jsZw0pzdUQsb6FlQiNnn4Oq1cTwy2EwIMTPNyOADQBsadbLpQ9M/hGCM3HcmgxNHEev0EmPkPh97qio5XrR9XIczgMNWOYc7iunrHf6Y2gGzx8AY52ARpcH2aWwfuLtNgkWFiexE7/iaHYxAABXw01dqEpk19MNMjvboz0pFTrjwFxUm+l9EUu280X/OmIM5yvWZ3VHsltbxH1AplMeDGqVz6EjsQi/RQhOXW6FNaWcRQsxcMzb4GESRyIKMMsmDtnI/E9VfDszn2paPjLmfaTj4XBZR2yQaycywmt4qReN7FxE8bhP4w8voaHTU/Q2fj+61Z9FTZRAtUQwbyZ3WFBauN4h8cBlFB+J4H9uWcj8TbcXxvPomXo9nxH2c7tQ6ippFhUTPW4DtSz1roQ2FN2yj8qC7jBBC5JIEHycpRp9AQfV3Cjx9l+HikQxlSOJdHlq6Rr6RAonyGBrDazOBH+fjZvs1vX5QtkkwWIDjaHR8+BSD8npRwQJ6Kwx6q+H0D2ikmM4KaqwgeMMmPkOhU50GXWvs1rYU33QapVCuUwGc8vuJlLmJ+8f34Q78kzo5UdjwO47jYDoxr4jMcCqEmOok+Bgn5UpkRrStwB5wpzB14rV+2lYQz1Fbfn9ZlKEzGngMFClTNL9zEYETNgWP7xjSBGMUFNDzV8voLTeJlA6N3vTAbEcetkSY5eVEl9clRqQAjtvAHiP+O8ndpwd1ZjWsEQIGxahzfCTLEnfw72tBd/VMuLOvEELkAwk+xutkZsQZPCxSGRrDPHVTcfQorTNZ+LGqzMwEHspRiQXYVKK/zEm2X9NdB7bHRTDgh0j0VPOIYaICAbpmu+ir1Iya8cinwMM4NQJFFfjpqXIng4/xUjoRbJjRsdM1WqnxXb+j0e2d2B0dsqKtEGJKk+AjRcp00EolbvBqQKfV5OsaPMM06muV8xE1kxE8pCl/pZuuuQFazzSHNFHFiqDlqjPwt9kUPJ0YkdNz6WL6ykxiRUOPl68ZD7OyAru+MjmfR9TvwoxpzPipfdRITR568Mdy+hDeIVQKgQfgeAysJXWYvdWw9+CEJjYTQoh8MGODD+WoRBbDSDEaUP3ZhRGaNpQC5RpmJIej86LT6kT52x30C69S4FpB27KCIZdh+zXh+RAtdlGwxYdSio75LqKlQzMeg+TbxxEqoK/Gjx4wmZhhQ/Ia+tdYUcMkNE6uwTJeqQQekBgiHi3zYoTc+A97QYIPIcQUNSODD3XsBPW/C9Jb6aH5IhsVSG2uCgD0aZlvxejNHSdH1Aw5jkoEJU6+3YUnxiqAtjWJGVGtQI4LMw6u2TXYVSXJ57GQd1z9LyYlxYyHEEJMNzMy+LBbWvE81op/4Txazq1AT+AmqU/76a/Qo95MRsyIaI2283uWy2RTwjhuyrZP07Ho5H7538bklBTSUx/M6jkl8BBCzHQzMvhIamlnzp/L6Z3lpeV8GxWcQAZkAN3fITPVH87K5YA5SgbENbRvyfAHSn9n0+BhKHojhuM2iPzN+fTOMrPWchQ4rijeG8N3vBvHnsTkGEphLl6AVTp0AbpIMAv/BE4fzTLBz09ZGv/RLozuCM4E15MRQoh8MKODD7u9Hd+jWwjWzqF9aS3OJIKPk5mQsTIgpzuZERktA6IMjTFck00WhI5auP+8Ffuyczny9pMzpWYnoxFosvH8aSvOJEd2KNOkr76YrmGmdHf1atx9mb2edGU6jLiDevMY1iTXmBFCiFxLeUKKo0eP8qEPfYiysjL8fj9nnXUWL774YvJ1rTW333471dXV+P1+1qxZw969e9Na6HynT+8PMg6JIMQB99CHMnMTeAB0zXYRv2IlnQ0jT/ueboHjiqpNNoV7uyY3pNQwMc+Yjz53CbHCzDdtaaXQxtCHNLEkSN0hhDgppeCjvb2diy++GLfbzR/+8Ad27tzJv/7rv1JScqrD3j333MO3v/1t7rvvPjZv3kxBQQFr164lMsG1OaYarRXaURNKDhgujelxhjxUDtcf65kDx97ipqshe+cMHbbw/W4Lzksjr/A7HsrtoueMMlqXB4kWZSECUCM8hNQdQohBUmp2+drXvkZtbS0PPPBAcltDw6m7ktaab33rW3zxi1/k6quvBuAnP/kJlZWVPPzww7z//e9PU7HTS3f3ULFV0zvLR+cSa9J9P6YTrUafdX6ilKMoOKrxtQ/N6hQc7GYiuR7l9mA01KL9iSlItcsgXmAMGwAYcVA2417Nd/CJQBsMO9w2E5Sl8Z3ow+iJomPjXAUwz0zXukMIMTEp/ab+7W9/y8qVK3nPe95DRUUF55xzDj/4wQ+Srx84cIDGxkbWrFmT3FZUVMSqVavYuHHjsMeMRqOEw+FBj2yz29sJ/XIzs3+9H3dHfo88ySalMxN4QOLGXfZqL8H/3jLkMdGMh+H3EV5ezonzizlxfjEt5xQSCw2fejBjGk+PgznWRGBDCs6AJpUJFTNlroiN2n0Q+/W9E17dN9ema90hhJiYlKrPN954g+9///ssXLiQxx57jE984hN8+tOf5sc//jEAjY2NAFRWVg56X2VlZfK10919990UFRUlH7W1tRO5jsnTGt3bR+lOTcErPuia/GqwE+n7kQuJLASUvq6TD19bZtoLTp6rZLeD60T41Ic08JHqMb1eXPPm4syvxfKq6dns4eSu3086TOu6QwiRspSaXRzHYeXKlXz1q18F4JxzzmHHjh3cd999XH/99RMqwK233sr69euTz8PhcM4qETscpvgXL1JaUsKBTy4kGoqP/aYRaJ3o96FUaqNfckE5MGtbN2rb68lt7r85l+aS9GeBlA2zXgzDy7uxJzN8dgAjWEDHyiqiIencma+me90hhEhNSsFHdXU1S5cuHbRtyZIl/OY3vwGgqqoKgKamJqqrq5P7NDU1cfbZZw97TK/Xi9ebvZEUY9GWhe7pofANTXfMR98cGyYRhMCpH/PZWml2LI6lcB/24u5WeNs1rj6N2d6DbQ3oAJHmjI1yFP4mjafLwWzrxrIm3q/G8PkwqirASCTunJA/sdJstjvmZrnvx1Q2E+oOIcT4pRR8XHzxxezevXvQtj179lBfXw8kOpBVVVWxYcOGZIURDofZvHkzn/jEJ9JT4ixwensp+dU2ykJBDv7DIvqWTK8MiI6YzH20B+OVfeA4aK2x45ntZGvEoPK5dvTuA1iT7DRplJfRdlENtifxgWoFOttddRRoTmZahl/rRZwyU+oOIcT4pBR83HLLLVx00UV89atf5b3vfS9btmzh/vvv5/777wdAKcXNN9/MXXfdxcKFC2loaOC2226jpqaGa665JhPlzxgdjeIAwSMax+sjUjH5UTB5kwHRCiNi4fSOPEumtzVGwRE/sSJFPDTxNIiyFd428HRpjM4erGh0yD5mSQkqOHT20ZHY5UXYHoUz+W45k6NIZIhOy4A4ruE7oxq2RqXQ0qRsjbvLwuyOQZqaqHJlJtUdQoixpRR8nH/++Tz00EPceuutfPnLX6ahoYFvfetbXHfddcl9Pve5z9HT08ONN95IR0cHb3nLW/jjH/+Iz+dLe+EzTUejlP7mZcoKCjh8/UJ6lk9yBtQ8y4CMxty8k4rtbnrftoymCyaeVjCjULOhBX3wKFbf8CM1rCV1dC4Y/wI72gAnX+bmPS0DAmD5Fbb7tLV/NLgiGlff+FMkrl4b98v7cXp7cSbRTJUPZlrdIYQYndI6v8ZjhMNhioqKuIyrcalc/7RNUC4XHe9fSdtSRbx4cmvADAw+sp0BcSyF0eLB02Ew91fN2Lv3jfkefdEKTpx7Kith+RmUCTFiCk/nyENyXb2aiscOYh09NvgFpTDPmE+8IkTbEh89sxWeMHjCmf1zVDagwd2nMWNpPJcG1f9PKVZg4LiHfrmuyAjn1GBYQ7d7OuK4Xng9J8NrLR3nKR6hs7OTwsLCrJ9/IvKx7hBiJkml3siX3495TVsWpb/dSdmGAEfeO4+u89KTAdFkNwDRPS7m/Xcfrn3HcDo6x/Uec9tuqnefykpEVzRw7BIvun+hO/8JTdWjB9Ej3SBtGyvcPWSz8njY/bFy3nXZZopcfXiNON/b8jbKNmbwpqEHBADpjnEGZkBG+E4tr0r2UxnIjGmUI/1GhBAzhwQf42SHw9DVReDEXHqavNghGxWYXBCilM5oPxDDcPD64hj9aYlwzMAVjmCfODHuYziRCAwILDxtVXjavWAkCuztcHBa21L+da6Uwim2uLzw1IRihaU99FaUDNnXFUkxI+IkMhvKSdzwk51RdaJfRsZu8iczWv3n0KcHIophVwQ+fZtyNGbEwYxY5FliUggh0kKCj1RoTcljeyh9roBjV9XSeeEkO6A6iQCk/9Cg0huEFASiXL9gE7NcXQA81bGYg4ULJtflZNcbzD5eiOovqI7GsIfpRDoR6xY9zdGGocHHb/afjfOnQoxx9rk041CytQU6wkTPrCVaMuDPPAv3cjOayK5YPoXjSv3Tdoct3DsPovsiOGn6bIUQIp9I8JEiu7UNWtsItMyms92D9tko/8RHIugBP3vVJOcyNwwHt8dKBjBF/ghzPS0EVJQXe+dxuKcYw3Imdf89PROSTnXuVurcrUO2by2t40BJUTL4UBa4u/XIU79rUN292G0dGNEaVJZ7pyrNoCzLkAzIwHIOmMJeORojrjH74jidYfQU72QqhBAjkeBjgoqe3E/Ry0U0XTaL9ovTNAxyrNnFx8iMlBb1cF39C3iNxLwkPhWj2OjlT+EzeeIbF1G8pwe1681s/PhPq3dWbufVd59qKnq+sYHYY7Nw9eX3lZjRRDCRGBY8TF+PuMYVOdXfw9MRx7X7cCLjIYGHEGIak+BjguwTJ+DECQJLy2jrdoFbo7yTC0L0cB0CBlAkRsoYhoPpGtpxodjXxyLvMQqMwan6DitAyWthnJd2Zj/wMEzMslKU57SOpB435jg/r9nudmYXtSefd8V9bArOQpv9n5cDrr5RMiE5onRidI1zWrkSq+lqXFGNq/fU92hE4thtHeBM7Tk9hBBiLBJ8TFLhxjdZfKCMEyuLaHurndE5PLRW4Gjqq9q4qupVzNNCiZDZh8+Y3FTw6eaqnMXrX6ynev7gTq4K+HDN8KuVjuWtxXtwXWMT629O2d9ZTvufq/F25Fn0MYLQkSjePY1g22jbAd0fgMTi2BJ4CCFmAAk+JslqbILGJoKzz6e1zwVuB+XOzHAKw9QopZnl6+Z8/xsZOcdYlNuDMcxspNqycLq7B03jaoZCOBUlnLHkCDfX/TltZahydfLO0m3J58/6zuA3BdWYkf7IT+vkui+5plWi34cRPxUYucLRofOeCDHNKa8XIxQc+oJlYXeGp8YS4CJtJPhIk4Jth1nYVknb0gCtb4mhzPT/Q2qoauHSWXup87Sk/djj1ffXZ9P8932Y5uAAK7K/kDPu2Zccxuuqm8Pr/1hDcX0Hf1/5WkbLtMR/jPPX7qA7nlhk7GBnCfHXyjAamzN63rHYHoVjQuGBPtzHTjUb6c4uJL8hZpquvzsbz42NeE4btrZnbw1Lv3Qw8UNOzBgSfKSJdbwRdbyRosB5tK4y0DhpC0CUoVEK6graeVtw59hvOI2BxvG5MHy+xNDNFH5hKJcL5fEkn4frXXzt7N8M6VdyT9GVUFKI0dMDgF1eyIXn7uHvK59PubypKjO7+eCszcnnm4vm82TJxQQK/MS95pCpzgcy7DRP7jVgplP6Jx1zdfRhvXkojScRIv8l647+XvLhBpONS35J0Bg8Xf7NxSvZUzoPo2voZIQAaI0TiUpfqGlGgo808+06zvwHq+iY56X1ovikm2CU0iyoPsF5pYc4w984oWOsLtzHX74wj64jK1j8ndZxTat+Uud7V9J85alAY+Hsg/jU0H4layt3cv9dFxOPLgLAF4hxTcke3CoxasPWBk6W1rxf4GvikY+HOfyBeQRDvfg9I/eDad0+i4qt6Ys+lKMpeLMLo6MbTBNtGnBi6PBhIaa71r8/H/uatuTzK2Zvwq88Q/b7P6XP89lvzSEcmTPscbp7fMz9DqjnX85YWUX2SfCRZtbRY5hHj1ESP5vWC9xgKjAmlgE5uQ7M/MIWrip6acJlmus+wfeW/4w/NZzF87++AHO/a+gcEoaJMoZmCNoXGdx/0U/GPMcK/0HuPe/giK+bygFNVgKQMrObry59eFz7fqLpQ7A1ff8MlANGS6f06RDTjnIl/p2Mu+5Yptm38penbR367/88r4enznx4xPO+Eovw8YdvptDlQtu29A2ZJiT4yBDPgWbqH55N92wPrRfGUZ7Ufl0rpVlQc4JFRU2sKDicljIt8R3jwU/YOH9zPgt/3I6zY1fiXG4PTf+wko6zhs4tsWTxm2k5N2Q3ABFCpI+rqpI9N8/DKrZZfF83zkuJ5l/l9XL0pvPoO6dvyHuuXbopLeeuNR1CNx5h99qzWfCAjfGX7Wk5rsgtCT4yxDp6DO/RY/hWnknrygCkGHygYHnJUa4p3pq2MtW6W7n//J/ywpkN/O+Gt+PZ0X8qj5uui/v4zwvHznBMVl4GICP9kBrvsGn5ISamOV1SyAevfIarCl/iH/+wDv9Lie3K48F3WQuvnvurjJ27xAzw2JJHObSwm2uf/SdK/5KxU4kskuAjw8yjLdT+sY6eSi9t51ko3+idppSR6ONRF2xnWeBo2sphD7jZ17g7OHSdhfH21QA4Ls1l83ek7VxjyacAxGh3E9o1tE+GHfLRV+NHD5NOHvT+uMZ/pAcjEut/ozPuFYOFyHeu2jns+1gtdkOE/xd6lbmuGN0f7eTwlRckdjA1n5/3v1kpS6nhwv++Rvaeu4p5v4ljPrlt7DeJvCXBR4ZZxxvxP9xIcNki2s8shnEEHxeWH+CK0KtpLYejDez+n/JVro5x9ePIpHwJQDxhA/v1vUO2u2rnQLV/zPcrS2McOp5Y80eIacauLOaz1z7CjUXHSPTXKGDbkH4c2RE0fDxz1kN0Luvj0j2fpfLJnBRDpIkEH9nS0k71X0qwfIm5KOIBaF9uo4KJfhaGqamvbKU6EKbBO/4l76eyXAUgNorbdl5N92ul1GwZfiSM7uoiuK/z1BTuI1BxG93Tm4liCiGG4VMufH/dzIGa1dT+KYbrifQ1TYvskeAjS+ymZoK/OjXplWvObLrm1mGfDD4MhysqX+ctBbtzVcScyFUAEv9LGQ1fH3kOErujE6T5RIi841VuNp3930RXxDmn6zPUPpHrEomJkOAjR3R3N7O2a6IHEpkQK6DY1VCV9uDDzoN+FWPJVADSYQe485W/IXo8MPgFrajbkV9r4AiRT1y1czhybR1dDQ7LvEcYbohsrrkwKby4mcP/7yKqN0YlAzLFSPCRI3ZHJ6FfbSbU/9xVWcFra6qgIr3nGdjXI59lIgA5YRVS/mCAwMMvDn1RZ2b9HSGmg+j8Cm7/5E/5u4J23MrMdXGGZSqD51b8CmuFzVmhT9MgGZApRYKPXBowWY7u7aP7+flc3/kRrlv2wqQ7nE6FjMfp0hWANFsh/nXXXxE+HmLR4V60TMssxLi46mtpumIOHWfAfPcJ3Mqb6yKNylQGJgaly0/QdNNFzHqlD+NpmQdkKpDgI0/Y4TB197yIUVrMr757LlecN7ngw9GJG/hUyHoMlI4A5EC0gpL/CFL9xMvoWCyNpRNieus+q5qvfv4/ucjXNWQNlnz29IoHiSy3OOehm1n4dK5LI8ZDgo88ouMxdGcY9fICPu76EO+a93LKGZCBGY+pFnicNFIAsrFnAf+972xse/TAJBr2svh4T2IRPSHEmFxz62i9uIbW5Yr57naCRjDXRUqJV7nxKjfVC0/Q/uHVFO3rw3juZZmKPY9J8JFnnEiE+m++hAoF+ek3LuCKi8cffNgYUzbjcbrhApD/3nc2c9eH0Z3h0d/saOzungyXUIjpo211DV/58g84yxOmzAiM/YY89b9n/ozOL9tc+qebWbTJHLoOjcgbEnzkIae3F2XbOPG6Cb1/qgceJz3XcwaPHV2CoxPXY+0qxGl5E6dHAgsh0slxKea6Oqkwp1bG43RFhp8iAwyv9PPKdxJ8iLz1ky0XsfTOY9D/66U82oItgYcQQkx5EnzkK9vGe8TN1w9fmdy0uvSNGTEJ2caehWxqayBwwI11rBFktIoQaeFqqCcyrxx9WnK0c4HCOw0Sps9E4I+dyzEP+9CO9PfIZxJ85CltWcy7dz+xnxQkt/3HZ97OW67aPahvx+mmQ5PLf/3lUhb/WzP1nXuxJfAQIm0Ovm82X/3ojxJ9qgaoMsNUm1O3r8dJNzx3A4vu6mJh5xtYUnfkNQk+8pjd1AxNp54XHKrmJy0XsyjQxAr/wdwVLEO29c1lV3c1gUMm9r4DuS6OENNOPKS5KtCNqU7/8eLJSXnSZUs0zra+ubgP+LB3yzwfU4EEH1NI3f9/P4f+NJcNHzqb773zP4f8epnq7vvTX3HGjzqoP/EG0kddCDFe7/vzJ1l0Xy/zmw5J3TFFpDSTk23b3HbbbTQ0NOD3+5k/fz7//M//jB44U6fW3H777VRXV+P3+1mzZg179w5dslykzmpswnn5dYKHFI91nsWfwonHy331uS7ahLwaqU1ewx87lxM8bOC8sgvreGOuiybSTOqO/OAOKx7uKebhniAP9wTZGo1hT8GlBl6KRpPX8HBPkMABN3rra1hHjua6aGKcUsp8fO1rX+P73/8+P/7xj1m2bBkvvvgiN9xwA0VFRXz6058G4J577uHb3/42P/7xj2loaOC2225j7dq17Ny5E59v6syYl89mP3SQV7ctTz4/8vYA//x/fkqxOXVGgnQ5fu7773dQ96dTy9HPOSS/WqYrqTvyQ/0vj/D9Le9OPj/4115evO6bFCl/DkuVmm4nwrt/czNzf39q9uK5bx6VumOKSSn4eP7557n66qu56qqrAJg7dy4PPvggW7ZsARK/XL71rW/xxS9+kauvvhqAn/zkJ1RWVvLwww/z/ve/P83Fn5msI0dRAyL8ovoL2dozlxJ3IvgIGRHmeZoz0izTagc5GCuf9HE6rQBF+zXquZeS26TymL6k7sgP1puHMN88lHxe3LCaZyMlzDK7AAipOGe4fcP0CZm8FruHN6zJ9y3psEMU7VGYT25LbpO6Y+pJKfi46KKLuP/++9mzZw9nnHEGL7/8Ms8++yzf/OY3AThw4ACNjY2sWbMm+Z6ioiJWrVrFxo0bpQLJkNKnD7Gp+Xy0mRjp0nqmh09+9BEWetPffPGVne+g5AdBlDW5YWxKQ/lOyXTMFFJ35KeKxw9zT/Pfo/sXrm1Z7uI3H/0GSzzpH/ny4f3vpuM7dRiTrTtsqHr1iNQdU1xKwccXvvAFwuEwixcvxjRNbNvmK1/5Ctdddx0AjY2Jm11lZeWg91VWViZfO100GiU6YA2OcHiMqbPFENbRY7iPHks+n2Wdx87eGtwq/UPNuo4WUv2n7ej45Bdsk8pj5pC6Iz9Zh4/gO3wk+bwifj4vf2g2HnU47ed6/UgVCx99CZ2GNZek7pj6Ugo+fvWrX/Gzn/2Mn//85yxbtoyXXnqJm2++mZqaGq6//voJFeDuu+/mzjvvnNB7xfC8Ow7z0j+fw1bPuWk/9sLDfWgrnvbjiulN6o6poWD7Yb5z2/uwMjDj2Lw3I7LKtEhSWo9/2b/a2lq+8IUvsG7duuS2u+66i5/+9Kfs2rWLN954g/nz57N9+3bOPvvs5D6XXnopZ599Nv/+7/8+5JjD/Xqpra3lMq7GpdwTvCwhxGRYOs5TPEJnZyeFhYWTPp7UHUJMf6nUGyn1Kurt7cUwBr/FNE0cJ9GxsaGhgaqqKjZs2JB8PRwOs3nzZlavXj3sMb1eL4WFhYMeQojpReoOIcRAKTW7/O3f/i1f+cpXqKurY9myZWzfvp1vfvObfOQjHwFAKcXNN9/MXXfdxcKFC5PD5WpqarjmmmsyUX4hxBQgdYcQYqCUgo/vfOc73HbbbXzyk5+kubmZmpoaPvaxj3H77bcn9/nc5z5HT08PN954Ix0dHbzlLW/hj3/8o4zTF2IGk7pDCDFQSn0+siEcDlNUVCTttkLkULr7fGSD1B1C5FbG+nwIIYQQQkyWBB9CCCGEyCoJPoQQQgiRVRJ8CCGEECKrJPgQQgghRFZJ8CGEEEKIrJLgQwghhBBZJcGHEEIIIbJKgg8hhBBCZJUEH0IIIYTIKgk+hBBCCJFVEnwIIYQQIqsk+BBCCCFEVknwIYQQQoiscuW6AKfTWgNgEQed48IIMUNZxIFT/x6nAqk7hMitVOqNvAs+urq6AHiW/81xSYQQXV1dFBUV5boY4yJ1hxD5YTz1htJ59tPGcRx2797N0qVLOXz4MIWFhbku0qSFw2Fqa2vlevLUdLqedF2L1pquri5qamowjKnROit1R36bTtcCcj3DSaXeyLvMh2EYzJ49G4DCwsJp8aWeJNeT36bT9aTjWqZKxuMkqTumhul0LSDXc7rx1htT4yeNEEIIIaYNCT6EEEIIkVV5GXx4vV7uuOMOvF5vrouSFnI9+W06Xc90upaJmG7XP52uZzpdC8j1TFbedTgVQgghxPSWl5kPIYQQQkxfEnwIIYQQIqsk+BBCCCFEVknwIYQQQoisysvg495772Xu3Ln4fD5WrVrFli1bcl2kMd19992cf/75hEIhKioquOaaa9i9e/egfS677DKUUoMeH//4x3NU4tF96UtfGlLWxYsXJ1+PRCKsW7eOsrIygsEg1157LU1NTTks8ejmzp075HqUUqxbtw7I/+/mmWee4W//9m+pqalBKcXDDz886HWtNbfffjvV1dX4/X7WrFnD3r17B+3T1tbGddddR2FhIcXFxXz0ox+lu7s7i1eRWVOx3gCpO/K57pB6I3P1Rt4FH7/85S9Zv349d9xxB9u2bWPFihWsXbuW5ubmXBdtVE8//TTr1q1j06ZNPP7448Tjca644gp6enoG7fcP//APHD9+PPm45557clTisS1btmxQWZ999tnka7fccgu/+93v+PWvf83TTz/NsWPHeNe73pXD0o7uhRdeGHQtjz/+OADvec97kvvk83fT09PDihUruPfee4d9/Z577uHb3/429913H5s3b6agoIC1a9cSiUSS+1x33XW89tprPP744zz66KM888wz3Hjjjdm6hIyaqvUGSN2Rz3WH1BsZrDd0nrngggv0unXrks9t29Y1NTX67rvvzmGpUtfc3KwB/fTTTye3XXrppfozn/lM7gqVgjvuuEOvWLFi2Nc6Ojq02+3Wv/71r5PbXn/9dQ3ojRs3ZqmEk/OZz3xGz58/XzuOo7WeWt8NoB966KHkc8dxdFVVlf7617+e3NbR0aG9Xq9+8MEHtdZa79y5UwP6hRdeSO7zhz/8QSul9NGjR7NW9kyZLvWG1lJ35DOpN9JXb+RV5iMWi7F161bWrFmT3GYYBmvWrGHjxo05LFnqOjs7ASgtLR20/Wc/+xnl5eWceeaZ3HrrrfT29uaieOOyd+9eampqmDdvHtdddx2HDh0CYOvWrcTj8UHf0+LFi6mrq5sS31MsFuOnP/0pH/nIR1BKJbdPpe9moAMHDtDY2Djo+ygqKmLVqlXJ72Pjxo0UFxezcuXK5D5r1qzBMAw2b96c9TKn03SqN0Dqjnwl9UZCuuqNvFpYrqWlBdu2qaysHLS9srKSXbt25ahUqXMch5tvvpmLL76YM888M7n9gx/8IPX19dTU1PDKK6/w+c9/nt27d/M///M/OSzt8FatWsWPfvQjFi1axPHjx7nzzjt561vfyo4dO2hsbMTj8VBcXDzoPZWVlTQ2NuamwCl4+OGH6ejo4MMf/nBy21T6bk538jMf7t/NydcaGxupqKgY9LrL5aK0tHRKfGejmS71Bkjdkc+k3khIV72RV8HHdLFu3Tp27NgxqJ0TGNROdtZZZ1FdXc3ll1/O/v37mT9/fraLOaorr7wy+d/Lly9n1apV1NfX86tf/Qq/35/Dkk3eD3/4Q6688kpqamqS26bSdyOmL6k78pfUG+mVV80u5eXlmKY5pOdzU1MTVVVVOSpVaj71qU/x6KOP8uSTTzJnzpxR9121ahUA+/bty0bRJqW4uJgzzjiDffv2UVVVRSwWo6OjY9A+U+F7OnjwIH/+85/5v//3/46631T6bk5+5qP9u6mqqhrS+dKyLNra2vL+OxvLdKg3QOqOfP6upN44JV31Rl4FHx6Ph/POO48NGzYktzmOw4YNG1i9enUOSzY2rTWf+tSneOihh3jiiSdoaGgY8z0vvfQSANXV1Rku3eR1d3ezf/9+qqurOe+883C73YO+p927d3Po0KG8/54eeOABKioquOqqq0bdbyp9Nw0NDVRVVQ36PsLhMJs3b05+H6tXr6ajo4OtW7cm93niiSdwHCdZYU5VU7neAKk7pkLdIfVGBuqNSXVXzYBf/OIX2uv16h/96Ed6586d+sYbb9TFxcW6sbEx10Ub1Sc+8QldVFSkn3rqKX38+PHko7e3V2ut9b59+/SXv/xl/eKLL+oDBw7oRx55RM+bN09fcsklOS758D772c/qp556Sh84cEA/99xzes2aNbq8vFw3NzdrrbX++Mc/ruvq6vQTTzyhX3zxRb169Wq9evXqHJd6dLZt67q6Ov35z39+0Pap8N10dXXp7du36+3bt2tAf/Ob39Tbt2/XBw8e1Fpr/S//8i+6uLhYP/LII/qVV17RV199tW5oaNB9fX3JY/z1X/+1Puecc/TmzZv1s88+qxcuXKg/8IEP5OqS0mqq1htaS92R73WH1BuZqTfyLvjQWuvvfOc7uq6uTns8Hn3BBRfoTZs25bpIYwKGfTzwwANaa60PHTqkL7nkEl1aWqq9Xq9esGCB/qd/+ifd2dmZ24KP4H3ve5+urq7WHo9Hz549W7/vfe/T+/btS77e19enP/nJT+qSkhIdCAT0O9/5Tn38+PEclnhsjz32mAb07t27B22fCt/Nk08+Oezf1/XXX6+1Tgybu+2223RlZaX2er368ssvH3Kdra2t+gMf+IAOBoO6sLBQ33DDDbqrqysHV5MZU7He0FrqjnyvO6TeyEy9obTWenK5EyGEEEKI8curPh9CCCGEmP4k+BBCCCFEVknwIYQQQoiskuBDCCGEEFklwYcQQgghskqCDyGEEEJklQQfQgghhMgqCT6EEEIIkVUSfAghhBAiqyT4EEIIIURWSfAhhBBCiKyS4EMIIYQQWfX/AQ0w8TUo2hSZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "TO_SEARCH = 8201 # 8201: temporal_mid_l\n",
    "f, ax = plt.subplots(1, 2)\n",
    "ax[0].imshow(atlasarr)\n",
    "ax[1].imshow(np.ma.masked_not_equal(atlasarr, TO_SEARCH).filled(0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((128, 128, 50), (91, 109, 91))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "heatmap = np.load(\"pos_mean_pet.npy\")\n",
    "heatmap.shape, atlas.get_fdata().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fix_heat_dim(heat):\n",
    "    # Both array and heatmap are numpy arrays\n",
    "    heat2 = np.full((128,128,128), 0, dtype=np.float64)\n",
    "    heat2[:,:,39:89] = heat.copy()\n",
    "    return heat2\n",
    "def fix_atlas_dim(atlas):\n",
    "    # Both array and heatmap are numpy arrays\n",
    "    atlas2 = np.full((128,128,128), 0, dtype=np.float64)\n",
    "    atlas2[19:110,10:119,19:110] = atlas.copy()\n",
    "    return atlas2\n",
    "atlas_fixed = fix_atlas_dim(atlas.get_fdata())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "heat_mri_pos = fix_heat_dim(np.load(\"resized_pos_mri.npy\"))\n",
    "heat_mri_neg = fix_heat_dim(np.load(\"resized_neg_mri.npy\"))\n",
    "heat_pet_pos = fix_heat_dim(np.load(\"pos_mean_pet.npy\"))\n",
    "heat_pet_neg = fix_heat_dim(np.load(\"neg_mean_pet.npy\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(data=[--, --, 7, --, --],\n",
       "             mask=[ True,  True, False,  True,  True],\n",
       "       fill_value=999999)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([1,2,3,4,5])\n",
    "b = np.array([5,6,7,8,9])\n",
    "m = np.ma.masked_not_equal(a, 3)\n",
    "np.ma.masked_where(m.mask, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>intensity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Precentral_L</td>\n",
       "      <td>2001</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Precentral_R</td>\n",
       "      <td>2002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Frontal_Sup_2_L</td>\n",
       "      <td>2101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Frontal_Sup_2_R</td>\n",
       "      <td>2102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Frontal_Mid_2_L</td>\n",
       "      <td>2201</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              name  intensity\n",
       "1     Precentral_L       2001\n",
       "2     Precentral_R       2002\n",
       "3  Frontal_Sup_2_L       2101\n",
       "4  Frontal_Sup_2_R       2102\n",
       "5  Frontal_Mid_2_L       2201"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\"AAL2_Atlas_Labels.csv\", names=[\"name\", \"intensity\"], index_col=0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []\n",
    "for _, row in df.iterrows():\n",
    "    to_search = row['intensity']\n",
    "    mask = np.ma.masked_not_equal(atlas_fixed, to_search)\n",
    "    pet_neg = np.ma.masked_where(mask.mask, heat_pet_neg).filled(0)\n",
    "    pet_pos = np.ma.masked_where(mask.mask, heat_pet_pos).filled(0)\n",
    "    mri_neg = np.ma.masked_where(mask.mask, heat_mri_neg).filled(0)\n",
    "    mri_pos = np.ma.masked_where(mask.mask, heat_mri_pos).filled(0)\n",
    "    d = {\n",
    "        'part': row['name'],\n",
    "        'Negative PET mean': pet_neg.mean(),\n",
    "        'Positive PET mean': pet_pos.mean(),\n",
    "        'Negative MRI mean': mri_neg.mean(),\n",
    "        'Positive MRI mean': mri_pos.mean(),\n",
    "        'Negative PET count': np.count_nonzero(pet_neg),\n",
    "        'Positive PET count': np.count_nonzero(pet_pos),\n",
    "        'Negative MRI count': np.count_nonzero(mri_neg),\n",
    "        'Positive MRI count': np.count_nonzero(mri_pos),\n",
    "        'Negative PET sum': pet_neg.sum(),\n",
    "        'Positive PET sum': pet_pos.sum(),\n",
    "        'Negative MRI sum': mri_neg.sum(),\n",
    "        'Positive MRI sum': mri_pos.sum()\n",
    "    }\n",
    "    l += [d]\n",
    "results = pd.DataFrame(l)\n",
    "del l;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>part</th>\n",
       "      <th>Negative PET mean</th>\n",
       "      <th>Positive PET mean</th>\n",
       "      <th>Negative MRI mean</th>\n",
       "      <th>Positive MRI mean</th>\n",
       "      <th>Negative PET count</th>\n",
       "      <th>Positive PET count</th>\n",
       "      <th>Negative MRI count</th>\n",
       "      <th>Positive MRI count</th>\n",
       "      <th>Negative PET sum</th>\n",
       "      <th>Positive PET sum</th>\n",
       "      <th>Negative MRI sum</th>\n",
       "      <th>Positive MRI sum</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Precentral_L</td>\n",
       "      <td>0.000532</td>\n",
       "      <td>0.000576</td>\n",
       "      <td>0.000461</td>\n",
       "      <td>0.000220</td>\n",
       "      <td>3038</td>\n",
       "      <td>3038</td>\n",
       "      <td>3038</td>\n",
       "      <td>3038</td>\n",
       "      <td>1114.707249</td>\n",
       "      <td>1207.641298</td>\n",
       "      <td>966.749000</td>\n",
       "      <td>461.281862</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Precentral_R</td>\n",
       "      <td>0.000427</td>\n",
       "      <td>0.000417</td>\n",
       "      <td>0.000250</td>\n",
       "      <td>0.000168</td>\n",
       "      <td>2747</td>\n",
       "      <td>2747</td>\n",
       "      <td>2681</td>\n",
       "      <td>2681</td>\n",
       "      <td>894.622546</td>\n",
       "      <td>873.637698</td>\n",
       "      <td>524.047146</td>\n",
       "      <td>352.063251</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Frontal_Sup_2_L</td>\n",
       "      <td>0.000673</td>\n",
       "      <td>0.000678</td>\n",
       "      <td>0.000188</td>\n",
       "      <td>0.000074</td>\n",
       "      <td>4472</td>\n",
       "      <td>4472</td>\n",
       "      <td>2321</td>\n",
       "      <td>2321</td>\n",
       "      <td>1411.819159</td>\n",
       "      <td>1421.030193</td>\n",
       "      <td>393.742026</td>\n",
       "      <td>154.649817</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Frontal_Sup_2_R</td>\n",
       "      <td>0.000678</td>\n",
       "      <td>0.000649</td>\n",
       "      <td>0.000242</td>\n",
       "      <td>0.000086</td>\n",
       "      <td>4689</td>\n",
       "      <td>4689</td>\n",
       "      <td>2541</td>\n",
       "      <td>2541</td>\n",
       "      <td>1422.782633</td>\n",
       "      <td>1361.921047</td>\n",
       "      <td>506.841342</td>\n",
       "      <td>180.509450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Frontal_Mid_2_L</td>\n",
       "      <td>0.000763</td>\n",
       "      <td>0.000788</td>\n",
       "      <td>0.000188</td>\n",
       "      <td>0.000094</td>\n",
       "      <td>4509</td>\n",
       "      <td>4509</td>\n",
       "      <td>3160</td>\n",
       "      <td>3160</td>\n",
       "      <td>1599.345071</td>\n",
       "      <td>1653.073211</td>\n",
       "      <td>395.204205</td>\n",
       "      <td>198.130331</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              part  Negative PET mean  Positive PET mean  Negative MRI mean  \\\n",
       "0     Precentral_L           0.000532           0.000576           0.000461   \n",
       "1     Precentral_R           0.000427           0.000417           0.000250   \n",
       "2  Frontal_Sup_2_L           0.000673           0.000678           0.000188   \n",
       "3  Frontal_Sup_2_R           0.000678           0.000649           0.000242   \n",
       "4  Frontal_Mid_2_L           0.000763           0.000788           0.000188   \n",
       "\n",
       "   Positive MRI mean  Negative PET count  Positive PET count  \\\n",
       "0           0.000220                3038                3038   \n",
       "1           0.000168                2747                2747   \n",
       "2           0.000074                4472                4472   \n",
       "3           0.000086                4689                4689   \n",
       "4           0.000094                4509                4509   \n",
       "\n",
       "   Negative MRI count  Positive MRI count  Negative PET sum  Positive PET sum  \\\n",
       "0                3038                3038       1114.707249       1207.641298   \n",
       "1                2681                2681        894.622546        873.637698   \n",
       "2                2321                2321       1411.819159       1421.030193   \n",
       "3                2541                2541       1422.782633       1361.921047   \n",
       "4                3160                3160       1599.345071       1653.073211   \n",
       "\n",
       "   Negative MRI sum  Positive MRI sum  \n",
       "0        966.749000        461.281862  \n",
       "1        524.047146        352.063251  \n",
       "2        393.742026        154.649817  \n",
       "3        506.841342        180.509450  \n",
       "4        395.204205        198.130331  "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.to_csv(\"heatmap_importance.csv\")\n",
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best for MRI:\n",
      "\tNeg: Cerebelum_10_L, Cerebelum_10_R, Cerebelum_7b_L, Cerebelum_7b_R, Cerebelum_9_R\n",
      "\tPos: Cerebelum_10_L, Cerebelum_7b_L, Cerebelum_7b_R, Cerebelum_10_R, Cerebelum_8_L\n",
      "Best for PET:\n",
      "\tNeg: Cerebelum_10_R, Cerebelum_10_L, Cerebelum_7b_R, Cerebelum_7b_L, Cerebelum_9_R\n",
      "\tPos: Cerebelum_10_L, Cerebelum_7b_R, Cerebelum_10_R, Cerebelum_7b_L, Cerebelum_9_R\n"
     ]
    }
   ],
   "source": [
    "pos_mri_best = results.sort_values(by=\"Positive MRI mean\").head()[\"part\"]\n",
    "neg_mri_best = results.sort_values(by=\"Negative MRI mean\").head()[\"part\"]\n",
    "pos_pet_best = results.sort_values(by=\"Positive PET mean\").head()[\"part\"]\n",
    "neg_pet_best = results.sort_values(by=\"Negative PET mean\").head()[\"part\"]\n",
    "\n",
    "print(f\"Best for MRI:\\n\\tNeg: {', '.join(neg_mri_best)}\\n\\tPos: {', '.join(pos_mri_best)}\")\n",
    "print(f\"Best for PET:\\n\\tNeg: {', '.join(neg_pet_best)}\\n\\tPos: {', '.join(pos_pet_best)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.000532\n",
       "1    0.000427\n",
       "2    0.000673\n",
       "3    0.000678\n",
       "4    0.000763\n",
       "Name: Negative PET mean, dtype: float64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results['Negative PET mean'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'numpy.ma.core.MaskedArray'>\n",
      "[   0. 2001.]\n",
      "[   0. 2002.]\n",
      "[   0. 2101.]\n",
      "[   0. 2102.]\n",
      "[   0. 2201.]\n",
      "[   0. 2202.]\n",
      "[   0. 2301.]\n",
      "[   0. 2302.]\n",
      "[   0. 2311.]\n",
      "[   0. 2312.]\n",
      "[   0. 2321.]\n",
      "[   0. 2322.]\n",
      "[   0. 2331.]\n",
      "[   0. 2332.]\n",
      "[   0. 2401.]\n",
      "[   0. 2402.]\n",
      "[   0. 2501.]\n",
      "[   0. 2502.]\n",
      "[   0. 2601.]\n",
      "[   0. 2602.]\n",
      "[   0. 2611.]\n",
      "[   0. 2612.]\n",
      "[   0. 2701.]\n",
      "[   0. 2702.]\n",
      "[   0. 2801.]\n",
      "[   0. 2802.]\n",
      "[   0. 2811.]\n",
      "[   0. 2812.]\n",
      "[   0. 2821.]\n",
      "[   0. 2822.]\n",
      "[   0. 2831.]\n",
      "[   0. 2832.]\n",
      "[   0. 3001.]\n",
      "[   0. 3002.]\n",
      "[   0. 4001.]\n",
      "[   0. 4002.]\n",
      "[   0. 4011.]\n",
      "[   0. 4012.]\n",
      "[   0. 4021.]\n",
      "[   0. 4022.]\n",
      "[   0. 4101.]\n",
      "[   0. 4102.]\n",
      "[   0. 4111.]\n",
      "[   0. 4112.]\n",
      "[   0. 4201.]\n",
      "[   0. 4202.]\n",
      "[   0. 5001.]\n",
      "[   0. 5002.]\n",
      "[   0. 5011.]\n",
      "[   0. 5012.]\n",
      "[   0. 5021.]\n",
      "[   0. 5022.]\n",
      "[   0. 5101.]\n",
      "[   0. 5102.]\n",
      "[   0. 5201.]\n",
      "[   0. 5202.]\n",
      "[   0. 5301.]\n",
      "[   0. 5302.]\n",
      "[   0. 5401.]\n",
      "[   0. 5402.]\n",
      "[   0. 6001.]\n",
      "[   0. 6002.]\n",
      "[   0. 6101.]\n",
      "[   0. 6102.]\n",
      "[   0. 6201.]\n",
      "[   0. 6202.]\n",
      "[   0. 6211.]\n",
      "[   0. 6212.]\n",
      "[   0. 6221.]\n",
      "[   0. 6222.]\n",
      "[   0. 6301.]\n",
      "[   0. 6302.]\n",
      "[   0. 6401.]\n",
      "[   0. 6402.]\n",
      "[   0. 7001.]\n",
      "[   0. 7002.]\n",
      "[   0. 7011.]\n",
      "[   0. 7012.]\n",
      "[   0. 7021.]\n",
      "[   0. 7022.]\n",
      "[   0. 7101.]\n",
      "[   0. 7102.]\n",
      "[   0. 8101.]\n",
      "[   0. 8102.]\n",
      "[   0. 8111.]\n",
      "[   0. 8112.]\n",
      "[   0. 8121.]\n",
      "[   0. 8122.]\n",
      "[   0. 8201.]\n",
      "[   0. 8202.]\n",
      "[   0. 8211.]\n",
      "[   0. 8212.]\n",
      "[   0. 8301.]\n",
      "[   0. 8302.]\n",
      "[   0. 9001.]\n",
      "[   0. 9002.]\n",
      "[   0. 9011.]\n",
      "[   0. 9012.]\n",
      "[   0. 9021.]\n",
      "[   0. 9022.]\n",
      "[   0. 9031.]\n",
      "[   0. 9032.]\n",
      "[   0. 9041.]\n",
      "[   0. 9042.]\n",
      "[   0. 9051.]\n",
      "[   0. 9052.]\n",
      "[   0. 9061.]\n",
      "[   0. 9062.]\n",
      "[   0. 9071.]\n",
      "[   0. 9072.]\n",
      "[   0. 9081.]\n",
      "[   0. 9082.]\n",
      "[   0. 9100.]\n",
      "[   0. 9110.]\n",
      "[   0. 9120.]\n",
      "[   0. 9130.]\n",
      "[   0. 9140.]\n",
      "[   0. 9150.]\n",
      "[   0. 9160.]\n",
      "[   0. 9170.]\n"
     ]
    }
   ],
   "source": [
    "print(type(np.ma.masked_not_equal(atlas.get_fdata(), TO_SEARCH)))\n",
    "\n",
    "for p in df['intensity']:\n",
    "    print(np.unique(np.ma.masked_not_equal(atlas.get_fdata(), p).filled(0)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.11.0 64-bit",
   "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.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}