phantom-generation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import skimage\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyphantom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hot-rod (Derenzo) Phantom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Produce derenzo phantom, with dimension 200x200\n",
      "nlayers: [18 12  9  6  5  5]\n"
     ]
    }
   ],
   "source": [
    "phantom=pyphantom.phantoms.generator.get_phantom('derenzo',shape=(200,200))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4gAAAMzCAYAAADkvTx4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbzElEQVR4nO39f5ScZX0w/r9nFrPhR3ZJwPzYGiCiVSoJKOg2Yi2UVEgslUeqgvE0IAVtAwKplSc98kusQbQWtTxgezTYI/jr8wC2FOkBFFIfQ4TQPAGrKeGLkBYCT4HskND82Jn7+0d2t8zssrOT+XXPva/XOXNOZua+rut93zOz+PZ6X/eVS5IkCQAAACa9fLsDAAAAIB0kiAAAAESEBBEAAIAhEkQAAAAiQoIIAADAEAkiAAAAESFBBAAAYIgEEQAAgIiQIAIAADBEgggAAEBESBABAABSb82aNXHaaadFX19f5HK5uP3228veT5IkLr/88pgzZ07sv//+sWjRonjsscdqHkeCCAAAkHI7duyIY445Jq6//vox37/22mvjK1/5Stx4442xbt26OPDAA+OUU06JnTt31jROLkmSpBEBAwAA0Hy5XC5uu+22OP300yNi7+xhX19f/Omf/ml88pOfjIiIgYGBmDVrVtx0001x5plnTrjv/ZoRcLOVSqV4+umnY9q0aZHL5dodDgAA0AJJksRLL70UfX19kc9PrBhy586dsXv37iZHtm+SJBmVz3R3d0d3d3dN/TzxxBOxdevWWLRo0chrvb290d/fH2vXrs1+gvj000/H3Llz2x0GAADQBlu2bInXve51VY/buXNnHLr//rGjBTHti4MOOii2b99e9toVV1wRV155ZU39bN26NSIiZs2aVfb6rFmzRt6bqI5MEKdNmxYREb9x8HHRlevIUwAAAGpUTAbjX7etH8kHqtm9e3fsiIgLuvJR25xc8+2KiL/evj22bNkSPT09I6/XOnvYaB2ZXQ1Pw3bl9pMgAgDAJFPrMrPuiOhO29K0oVvB9PT0lCWI+2L27NkREfHss8/GnDlzRl5/9tln49hjj62pL3cxBQAAMm2/XJLKR6PMmzcvZs+eHffee+/Ia4VCIdatWxcLFy6s7Vo1LCoAAACaYvv27bF58+aR50888URs2LAhZsyYEYcddlhcfPHF8dnPfjbe+MY3xrx58+Kyyy6Lvr6+kTudTpQEEQAAIOUeeuihOOmkk0aer1ixIiIili1bFjfddFN86lOfih07dsT5558f27Zti3e9611x1113xdSpU2sapyP3QSwUCtHb2xvzp/dbgwgAAJNEMRmMR15cFwMDAxNatzecN6x8TS6mpmwN4s4kiVV7kgmfS6tYgwgAAEBESBABAAAYoj4TAADItP1yex9pktZEzAwiAAAAESFBBAAAYEhaZzYBAAAaoiu395EmXe0O4FWYQQQAACAiJIgAAAAMUWIKAABkmhLTiTODCAAAQERIEAEAABiixBQAAMi0rlySwhLTpN0hjMkMIgAAABEhQQQAAGCIElMAACDTunIR+6WuxDSdzCACAAAQERJEAAAAhigxBQAAMq0rH+m7i2k6b2JqBhEAAIC9JIgAAABEhBJTAAAg47pyKSwxbXcAr8IMIgAAABEhQQQAAGCIElMAACDTlJhOnBlEAAAAIkKCCAAAwBAlpgAAQKYpMZ04M4gAAABEhAQRAACAIUpMAQCATOvKp7DENGl3BGMzgwgAAEBESBABAAAYosQUAADItHykr8Q0r8QUAACANJMgAgAAEBFKTAEAgIzryqWvxDRt8QwzgwgAAEBESBABAAAYosQUAADItHx+7yNNUhbOiLTGBQAAQItJEAEAAIgIJaYAAEDGuYvpxJlBBAAAICIkiAAAAAxRYgoAAGSaEtOJM4MIAABAREgQAQAAGKLEFAAAyLSu/N5HmnS1O4BXkbLLBAAAQLtIEAEAAIgIJaYAAEDG5XN7H2mStniGmUEEAAAgIiSIAAAADFFiCgAAZFo+n4uulNV07p2pS9ocxWhmEAEAAIgICSIAAABDlJgCAACZ5i6mE2cGEQAAgIiQIAIAADBEiSkAAJBpXfm9jzTpancAryJllwkAAIB2kSACAAAQEUpMAQCAjHMX04kzgwgAAEBESBABAAAYosQUAADINHcxnbiUXSYAAADaRYIIAABARCgxBQAAMi6fy0U+ZbcNzScREUm7wxjFDCIAAAARIUEEAABgiBJTAAAg0/K59G1Mn7Z4hplBBAAAICIkiAAAAAxRYgoAAGRaV37vI0260ncD04gwgwgAAMAQCSIAAAARocQUAADIuFw+F/mU3TY0p8QUAACANJMgAgAAEBESRAAAAIZYgwgAAGRaPrf3kSZpi2eYGUQAAAAiYh8SxDVr1sRpp50WfX19kcvl4vbbby97P5fLjfn4whe+MHLMEUccMer9a665pu6TAQAAYN/VXGK6Y8eOOOaYY+KjH/1ovP/97x/1/jPPPFP2/Ic//GGce+65ccYZZ5S9/pnPfCbOO++8kefTpk2rNRQAAICq8vm9jzTJp3Sbi5oTxMWLF8fixYtf9f3Zs2eXPf/BD34QJ510Urz+9a8ve33atGmjjgUAAKB9mppHP/vss/GP//iPce65545675prrolDDjkk3vrWt8YXvvCFGBwcbGYoAAAAVNHUu5h+85vfjGnTpo0qRf3EJz4Rb3vb22LGjBnx05/+NFauXBnPPPNMfOlLXxqzn127dsWuXbtGnhcKhWaGDQAAZEg+n4t8ym4bmpkS01p84xvfiKVLl8bUqVPLXl+xYsXIvxcsWBBTpkyJj33sY7Fq1aro7u4e1c+qVaviqquuamaoAAAAk17TSkz/+Z//OTZt2hR/9Ed/VPXY/v7+GBwcjF/96ldjvr9y5coYGBgYeWzZsqXB0QIAANC0GcSvf/3rcdxxx8UxxxxT9dgNGzZEPp+PmTNnjvl+d3f3mDOLAAAA1biL6cTVnCBu3749Nm/ePPL8iSeeiA0bNsSMGTPisMMOi4i9awS///3vx1/+5V+Oar927dpYt25dnHTSSTFt2rRYu3ZtXHLJJfGRj3wkpk+fXsepAAAAUI+aE8SHHnooTjrppJHnw+sJly1bFjfddFNERHznO9+JJEnirLPOGtW+u7s7vvOd78SVV14Zu3btinnz5sUll1xSti4RAACA1sslSZLSyc1XVygUore3N+ZP74+uXFPvswOkyPov7yh7ftxFB7ak7WQdu1PjbufYnRp3O8duZ9xA5ykmg/HIi+tiYGAgenp6qh4/nDc8cOx+cVBXuu5iur2YxG9uGJzwubRKyipxAQAAaBcJIgAAABHR5H0QAQAA2i2fz0U+n64S08zcxRSgVSrXFl12zsby91cvKHteufbole3raVtr+0bGXa19I+OubO+a1T62a1b72K28ZgBUp8QUAACAiDCDCAAAZFw+F5FP2dRYvtTuCMaWsssEAABAu9gHEegYnbrXWqeO3alxt3PsTo27nWPbBxGoxb7ug/jQca+Jg/ZL101qtg8mcfz6PanbB1F2BQAAZFoq72Ka0lrOlIYFAABAq0kQAQAAiAhrEIEMScs6qE5Z+5WmsTs17naO3alxt3psIFv2dQ3iv/RPiWkpW4P40mASb123O3VrEM0gAgAAEBESRAAAAIaozwQAADItl89FLmV3Mc2ldKpOggh0rMq1SZeds7H8/dULRv5duW6plra1tq+nba3tGxl3ZXvXrPaxXbPax2503ADUJ6V5KwAAAK1mBhEAAMi0fD59G9OnLZ5htrmANqkso6pGGVV1abmNfydtAZCWsTs17naO3alxt3psIFv2dZuLjSd0p3KbiwX/Z5dtLgAAAEgn028AAECm5fO5yKfsLqZpLTFNaVgAAAC0mjWI0CLVbt1ezdXj3Bae6jp5DVW7xu7UuNs5dqfG3c6xrTEEarGvaxAf/a2pqVyDePQ/70zdGkTZFQAAkGm5fC5yKSsxzaW0ljOlYQEAANBqEkQAAAAiQokpkGGvXNtUueZz/SvWdEaMXvdUbc3oeO3raVtr+0bGXdneNat9bNes9rFrbQuwL/L59N01NG3xDEtpWAAAALSaBBEAAICIUGIKAABknLuYTpx9EKFF7IPYXp26R1w7x+7UuNs5dqfG3c6x7YO4V+V1qMdkvYZMDvu6D+IvTz4glfsgvvnel1O3D2JK81YAAABazfQbAACQaflcLvIpKzFNWTgjzCACAAAQEdYgQtvUut7EmpLmauT6rU5Z+5WmsTs17naO3alx70v7LKp3Xfp4rraXJBm2r2sQ/23RgTHtNemasntpTxK/fs+O1K1BlF0BAACZlsun766haYtnWErDAgAAoNUkiAAAAESEElNoqFrWFda7hqeRfU9G1db/rK+yhme89vW0rbV9I+OubN/MuCvbu2a1t8/qNcuqZq45rFTr5wWTQb4rF/mudK1BzJfaHcHYzCACAAAQERJEAAAAhtjmAupQT8lQtduQN7NvRrP9QHvH7tS42zl2p8a9L+07UStLSmvlvxF0sn3d5uL/995pqdzm4vX/+FLqtrkwgwgAAEBESBABAAAYoj4TAADItFxXLnIpu4tpLqV3MZUgAkT1NTjV1lCN176ettXad2rc7Ry7U+Nu59jNjBuAdFFiCgAAQESYQQQAALIun4tcPl0lpmmdqktpWAAAALSaGUSAMVTbw2z9OPuI1dO21vb1tK3WvplxV7Z3zWpv38nXDIDaFYvFuPLKK+Nb3/pWbN26Nfr6+uLss8+OT3/605HLNW52VIIIAABkWr4rF/mU3cU0X+NdTD//+c/HDTfcEN/85jfjLW95Szz00ENxzjnnRG9vb3ziE59oWFwSRAAAgJT76U9/Gu973/vive99b0REHHHEEfHtb387fvaznzV0HGsQAQAA2qRQKJQ9du3aNeZx73znO+Pee++Nf/u3f4uIiP/7f/9v/OQnP4nFixc3NJ5ckiRJQ3tsgUKhEL29vTF/en905UyCkh6V63LGU+uanGb2TXXV9nlrVtt2jt2pcbdz7E6Nu91jZ1G1dZrtdLU1onSwYjIYj7y4LgYGBqKnp6fq8cN5w1MfnB49r0lXiWlhTxKHfe/FUa9fccUVceWVV456vVQqxZ//+Z/HtddeG11dXVEsFuMv/uIvYuXKlQ2NS3YFAADQJlu2bClLdru7u8c87nvf+17cfPPNccstt8Rb3vKW2LBhQ1x88cXR19cXy5Yta1g8EkQAAIA26enpmdBs6J/92Z/F//yf/zPOPPPMiIiYP39+PPnkk7Fq1SoJIgAAwETlchG5fLpKTHO52lb6vfzyy5HPl99CpqurK0qlGm+HWoU1iDCOWtb9RdS3ZqdT+mZsr7zmnbT2Ky1jd2rc7Ry7k+KejKxBhObY1zWIWz40PXqmpOv+nIXdpZj73RcnfC5nn3123HPPPfG1r30t3vKWt8S//Mu/xPnnnx8f/ehH4/Of/3zD4pJdAQAApNxXv/rVuOyyy+JP/uRP4rnnnou+vr742Mc+FpdffnlDx5EgAgAAmZbrykWuK2UlpjXGM23atLjuuuviuuuua05AQ9I1zwoAAEDbmEGEV6h3zcj6V6zrqFzT0Sl9j9U/o413zatdz2qfVy3t62lbrX0z465s75rV3j7N14wxrmfFNWvlmkRrDoFaSBABAIBMy0KJaasoMQUAACAiJIgAAAAMUWIKsA/GW19UbX1PtbVJtbSvp22t7RsZd2V716z29mm+ZozWyjWJ1hzCaLmufOS60jU3lutqdwRjS9dVAgAAoG0kiAAAAESEElOAhhivhKtyu4HKY6uVf43Xvp621do3M+52jt2pcbdzbCWKjVet5LSRfQPuYloLM4gAAABEhAQRAACAIUpMAQCATFNiOnESRHiFeteE1HO7+7T2Te0q13JV3r6+2pYBtbSvp2219s2Mu7K9a1Z7+3ZeMxrPNQbSQokpAAAAEWEGEQAAyLhcPoUlpvl0xTPMDCIAAAAREZFLkiRpdxC1KhQK0dvbG/On90dXziQoE1e55qaaWteE1NL/ZOmb6vvTNbN9u9pO1rE7Ne5GtM+iWv+bMZ40Xc+snheTQzEZjEdeXBcDAwPR09NT9fjhvOHZj8+Onu50zY0VdpVi1o1bJ3wurSK7AgAAMi3XlY9cV7oSxFxXuyMYW7quEgAAAG0jQQQAACAirEFkEhhvj7Jqrq5zH7HJ2De169R1Z40cu1PjbufY1hg2Xj1/G6tp59/OrJ4Xk9O+rkH8f5/4tVSuQXztV/4jdWsQ03WVAAAAaBsJIgAAABHhLqYAAEDWdeX2PtIkbfEMkSACtFi19UDra1hDWk/bau3raVtr+0bGXa29azax9pNBM9fmVWrl9c/qeQGtocQUAACAiNiHBHHNmjVx2mmnRV9fX+Ryubj99tvL3j/77LMjl8uVPU499dSyY1544YVYunRp9PT0xMEHHxznnntubN++va4TAQAAGEuuK5fKRxrVnCDu2LEjjjnmmLj++utf9ZhTTz01nnnmmZHHt7/97bL3ly5dGj//+c/j7rvvjjvuuCPWrFkT559/fu3RAwAA0DA1r0FcvHhxLF68eNxjuru7Y/bs2WO+94tf/CLuuuuuePDBB+P444+PiIivfvWrsWTJkvjiF78YfX19tYYE0FFGrRWrcc3OK9+vp22t7RsZd2X7ZsZd+b5rNrH2WdTKtXnVNHLtXlbPC2iPpqxBvO+++2LmzJnxpje9Kf74j/84nn/++ZH31q5dGwcffPBIchgRsWjRosjn87Fu3box+9u1a1cUCoWyBwAAwETkuvKpfKRRw6M69dRT4+/+7u/i3nvvjc9//vNx//33x+LFi6NYLEZExNatW2PmzJllbfbbb7+YMWNGbN26dcw+V61aFb29vSOPuXPnNjpsAACASa/h21yceeaZI/+eP39+LFiwII488si477774uSTT96nPleuXBkrVqwYeV4oFCSJAAAADdb0fRBf//rXx6GHHhqbN2+Ok08+OWbPnh3PPfdc2TGDg4PxwgsvvOq6xe7u7uju7m52qGTUeGuPamk7kfdr6T+rfVO7ate0cn3RK4+vp2219vW0TfPYnRp3s8cGyLR8Ln0b0+dTFs+Qphe+/vu//3s8//zzMWfOnIiIWLhwYWzbti3Wr18/csyPfvSjKJVK0d/f3+xwAAAAeBU1zyBu3749Nm/ePPL8iSeeiA0bNsSMGTNixowZcdVVV8UZZ5wRs2fPjscffzw+9alPxRve8IY45ZRTIiLiqKOOilNPPTXOO++8uPHGG2PPnj1xwQUXxJlnnukOpgAAAG1Uc4L40EMPxUknnTTyfHht4LJly+KGG26IjRs3xje/+c3Ytm1b9PX1xXve8564+uqry0pEb7755rjgggvi5JNPjnw+H2eccUZ85StfacDpMBlVllWNp9YSq1r6rrV/ffNqqt2yfrytDmppW619PW1rbd/IuKu1d80AJp80bkyftniG1ZwgnnjiiZEkyau+/0//9E9V+5gxY0bccssttQ4NAABAE6Vz8w0AAABarul3MQUAAGirrlxE2jam73r1qsx2kiDScaqtuxm3bZ3rgWrpv1P7ruy/mX0ztmrbkIx3DetpW/l+PW1rbd/IuNs5difFDQBjSVkaDQAAQLuYQQQAADLNXUwnzgwiAAAAEWEGESD1xltLVrlGtPLYauvQxmtfT9ta2zcy7mrtJ+s1A4CJkCACAADZ1pXb+0iTtMUzRIkpAAAAESFBBAAAYIgSU4AO88p1aJX7Uta712c9+23WM3Yj467WfrJcM0artndkrfu81uPqBn5+WT0vaKRcVz5yXemaG8t1Je0OYUzpukoAAAC0jQQRAACAiFBiCgAAZF0+hXcxzacsniG5JEnSWfw6jkKhEL29vTF/en905eS4k13lGp3x1LoWopa+a+1f3zRCtX31mtm+kWN3atztHNtvq/GqrRmtRzvX5mX1vJicislgPPLiuhgYGIienp6qxw/nDduufVP07N/VgggnrvBfxTj4U5smfC6tosQUAACAiFBiCgAAZF1XCktM0xbPEDOIAAAARIQ1iKRQJ6+fS8t6yMnSN9WlZb3cZBm7nXHTeLX+N2M8afoss3peTA77vAbxL9+czjWIf/rL1K1BlF0BAADZpsR0wpSYAgAAEBFmEEmBem+/vb7KLbZf2X8z+661f33X3jfVVbv+9fxe6v3sx2vfyLirtW9k3JXtW3nNaL6sXu+snhfQGBJEAAAg2/K59G1Mn7Z4higxBQAAICIkiAAAAAxRYgqQIaPWuNW4hq1sXWAdbWtt38i42zl2K+OerLK6RYPzqi5N50UH6srvfaRJVzp3G0zZVQIAAKBdJIgAAABEhBJTAAAg67py6duYPm3xDJEgAmRYtTU7leuDXnl8PW1rbV9P21rbNzLuyuObfc0mo3r3yh237zau+XRe+9C3NbrQEkpMAQAAiAgziAAAQNblU1himk9ZPEPMIAIAABAREbkkSdK5Acc4CoVC9Pb2xvzp/dGVMwmaNbXumVTLGoRm9l1r//quv29qV2190NWvWOMzao++GtrW2r6etrW2b2Tcle0bHfdk1Mw1bNU08/NwXo3n9zM5FZPBeOTFdTEwMBA9PT1Vjx/OG7b9zYLoOaCrBRFOXOHlYhx8/sYJn0uryK4AAIBscxfTCVNiCgAAQERIEAEAABhiDSIt0alr3Dp1PaS+mahq+/I1q21l+3p/m+2Ku9VjZ1E717BVU88aN+fVetYkTg77vAZx9bHpXIN4zobUrUE0gwgAAEBESBABAAAYoj4TAADItKQrF0nK7hqatniGSRBpinrWKKyvc0+ydvVd2X+n9l3ZfzP7ruy/0X1T3XjXrNp6uWrXu572zRw7zXEDQLspMQUAACAizCACAABZl8/tfaRJ2uIZYgYRAACAiDCDCMArjLcmtN41vPWsdW3k2I2Mu7J9rW0BIG0kiAAAQKa5i+nEKTEFAAAgIswgAvAKZWWfNZZHjiq9rKF9PW3TNLaSUgA6nQQRAADItnykr3YybfEMSWlYAAAAtJoEEQAAgIhQYkqTVFvTU0vbTum7Wv/6rv5+o/umPtWub+X2EJXH19K+nra1tm9k3JXH+04CpFPSlb67hiZd7Y5gbGYQAQAAiAgJIgAAAEOUmAIAANmWy0Xk01ViGrmUxTNEgsg+qVyDU00t63I6te9a+9d3uvtmtMrrfdk5G8vfr7IH4Hjt62lba/tGxl3Z3vesdtXWOVde72a6uoH7WDqv5mvkeQH/TYkpAAAAEWEGEQAAyLgkn4skZSWmaYtnmBlEAAAAIsIMIhNUbQ1O1fbjrBNoZt+V/Tez71r713e6+qa6amuP6tlntN49SuvZb7PZY1ObVq5xa+UaNudVP2sOoTUkiAAAQKYl+b2PNElbPMNSGhYAAACtJkEEAAAgIpSYArCPqq3/qVwzWnn8eO3raVutfTPjpvGqrXFrZN+t5Lzq7xtqocR04lIaFgAAAK0mQQQAACAilJgCAAAZl3TlIulK18b0aYtnmAQRgIaotk9lPfuh1rtH5nh7LtY7Nq2V1evvvIC0UGIKAABARJhBBAAAMs5dTCcupWEBAADQarkkSZJ2B1GrQqEQvb29MX96f3TlTIK2Q+WanWpqWYPQqX3X2r++09039au2n2Cz2rZ77Mmo1r+t40nT9XZe1TkvWq2YDMYjL66LgYGB6OnpqXr8cN6w9a53RM+B6cobCjsGY/apP5vwubRKuq4SAABAgykxnbiUhgUAAECrmUFkRKeWAqal71r713d2+mZs1a7hKz+TymNraVtr+3rasle1rUHq6ruN24o4r33o23lB5kgQAQCATFNiOnEpDQsAAIBWkyACAAAQEUpMJ7V6aver1eZPxr4r++/Uviv7b2bflf13Ut/UbrzvUr2/zXq+Cz7r6pq51qtSKz8f51U/50WnKOUiSimbGivl2h3B2FJ2mQAAAGgXCSIAAAARocQUAADIOHcxnTgJIgAtMWpd4DjramtpW2t7a4kmpp41vI3UyLVgrVzDVo3zqi6r5wVpl9K8FQAAgFarOUFcs2ZNnHbaadHX1xe5XC5uv/32kff27NkTl156acyfPz8OPPDA6Ovriz/8wz+Mp59+uqyPI444InK5XNnjmmuuqftkAAAAKpXy6XykUc1h7dixI4455pi4/vrrR7338ssvx8MPPxyXXXZZPPzww3HrrbfGpk2b4vd///dHHfuZz3wmnnnmmZHHhRdeuG9nAAAAQEPUvAZx8eLFsXjx4jHf6+3tjbvvvrvstb/+67+Od7zjHfHUU0/FYYcdNvL6tGnTYvbs2bUOD0BGjLdmp3LtUeWx1db7jNfeWiEAeHVNn9gcGBiIXC4XBx98cNnr11xzTRxyyCHx1re+Nb7whS/E4ODgq/axa9euKBQKZQ8AAICJGL6LadoeadTUu5ju3LkzLr300jjrrLOip6dn5PVPfOIT8ba3vS1mzJgRP/3pT2PlypXxzDPPxJe+9KUx+1m1alVcddVVzQwVAABg0mtagrhnz5744Ac/GEmSxA033FD23ooVK0b+vWDBgpgyZUp87GMfi1WrVkV3d/eovlauXFnWplAoxNy5c5sVOgAAwKTUlARxODl88skn40c/+lHZ7OFY+vv7Y3BwMH71q1/Fm970plHvd3d3j5k4Up9q+4rV0lbf1fvXd/X3O6lvGqvafmfV9hyrpb3PFmDySVJ419BJU2I6nBw+9thj8eMf/zgOOeSQqm02bNgQ+Xw+Zs6c2ehwAAAAmKCaE8Tt27fH5s2bR54/8cQTsWHDhpgxY0bMmTMn/uAP/iAefvjhuOOOO6JYLMbWrVsjImLGjBkxZcqUWLt2baxbty5OOumkmDZtWqxduzYuueSS+MhHPhLTp09v3JkBAABQk5oTxIceeihOOumkkefDawOXLVsWV155Zfz93/99REQce+yxZe1+/OMfx4knnhjd3d3xne98J6688srYtWtXzJs3Ly655JKyNYYAAACNkuSTSPJJu8Mok7Z4htWcIJ544omRJK9+MuO9FxHxtre9LR544IFah2UfVK7JqaaWdTn6rr9/fWen733pn3LV1rrWu37Y5wMAE5PSpZEAAAC0WlP3QQQAAGi3UgrvYpq2eIZJEDOk2m3eq7YfpySrU/uu7L+Zfdfav76z03dl/8oZ61ftGlZ+XpXH+wwAYN+kNG8FAACg1cwgAgAAmabEdOJSGhYAAACtZgYRgNSrtkbUthaNV7ZWvOL61rpGtx5XN/CzrbYdivOqn/OCzidBBAAAMq2UT6KUso3p0xbPMCWmAAAARIQEEQAAgCFKTAFIvWprkawHaq5WrgVr5Vov51U/50WnKOXSd9fQUq72Nv/xH/8Rl156afzwhz+Ml19+Od7whjfE6tWr4/jjj29YXBJEAACAlHvxxRfjhBNOiJNOOil++MMfxmtf+9p47LHHYvr06Q0dR4IIAACQcp///Odj7ty5sXr16pHX5s2b1/BxUjbRCgAA0FjDdzFN2yMiolAolD127do15jn8/d//fRx//PHxgQ98IGbOnBlvfetb42//9m8bfq1ySZKk8/6q4ygUCtHb2xvzp/dHV84k6Kup3Desmlrq6/Vdf//6zk7f+9I/ZEmtv5fxpOm35Lyqc160WjEZjEdeXBcDAwPR09NT9fjhvOHRDcfFtGnpyhteemkwjj52/ajXr7jiirjyyitHvT516tSIiFixYkV84AMfiAcffDAuuuiiuPHGG2PZsmUNiytdVwkAAGAS2bJlS1my293dPeZxpVIpjj/++Pjc5z4XERFvfetb49FHH5UgAgAA1KKUT+FdTIfi6enpmdBs6Jw5c+I3fuM3yl476qij4n//7//d0LhSdpkAAACodMIJJ8SmTZvKXvu3f/u3OPzwwxs6jhnEDpOWNVf6rr9/fWen71r7t7aFrMnqd9p5dZasnhcMu+SSS+Kd73xnfO5zn4sPfvCD8bOf/Sz+5m/+Jv7mb/6moeNIEAEAgEwr5fZtY/pmqjWet7/97XHbbbfFypUr4zOf+UzMmzcvrrvuuli6dGlD45IgAgAAdIDf+73fi9/7vd9r6hjWIAIAABARZhBTr3Jd02XnbJx429ULyp5X1ubru7F9V/bfqX1X9t/Mviv779S+a+2/Wt/QCFndE855Vee8mi+r55Vle+9imq7t39N2V9VhKQ0LAACAVpMgAgAAEBFKTAEAgIxL8ukr6UxSFs8wCSIAZEC964XH7buN62ad1z707bwaLqvnBWNJad4KAABAq5lBBAAAMq2US6KUS9ldTFMWzzAJIgB0oGaWvFWq7LuZJXHOq37Oq35ZPS+YCCWmAAAARIQZRAAAIOOK+b2PNElbPMNSGhYAAACtlkuSJJ2rI8dRKBSit7c35k/vj67c5JoErayJH0+tNev6rq/vWvvXd3b6rrV/60nYV6/8njVzTVStrq5jzVQr13rVynmN5rxar57zyqJiMhiPvLguBgYGoqenp+rxw3nDA//2tjhoWlcLIpy47S8V4zd//eEJn0urTK7sCgAAmHRK+b2PNElbPMNSGhYAAACtJkEEAAAgIpSYtl2nrrnSd/396zs7fdfaf7O/hwBAuVIuSd3G9GmLZ5gZRAAAACJCgggAAMAQJaYAAECmFfPp25g+bfEMkyC2WL377Kx/xV44leuSGtl3Zf+d2ndl/83su9b+9Z2dviv7b+VvEwCgkVKatwIAANBqZhABAIBMK0VEKdfuKMqV2h3AqzCDCAAAQERIEAEAABiixBQAAMi0UgrvYlpKWTzDUhoWAAAArSZBBAAAICKUmAIAABlXyqXwLqYpi2eYGUQAAAAiQoIIAADAECWmAABAphVzex9pkrZ4huWSJEnaHUStCoVC9Pb2xvzp/dGV6+wcd/2Xd9R0/HEXHajvlPZda//6zk7ftfbf7O8hk0Pl9+iycza2bOyrVy8oe97I76jzajznVbusnlcWFJPBeOTFdTEwMBA9PT1Vjx/OG+7a8tY4sKerBRFO3I5CMU6d+y8TPpdWUWIKAABARCgxBQAAMq6UT9/G9GmLZ5gSUwDIgGaWxLWz5M151c55NV5Wz6sT7WuJ6Z3/kc4S0yW/psQUAACAlDL9BgAAZJq7mE6cGUQAAAAiwhpEAMikWrdTGU+a1kQ5r+qcV/Nl9bw6wb6uQfyHp9O5BvG0vvStQZRdAQAAmVbK7X2kSdriGabEFAAAgIiQIAIAADBEiSkAZFBW1zU5r87ivEiLYi4XxVy6ajrTFs8wM4gAAABEhAQRAACAIUpMAQCATHMX04kzgwgAAEBESBABAAAYosQUAADItFIK72JaSlk8w8wgAgAAEBESRAAAAIYoMQUAADKtmNv7SJO0xTPMDCIAAAARIUEEAABgiBJTAAAg00q59N01tJSucEaYQQQAACAiJIgAAAAMUWIKAABkWjGXi2LKSkzTFs8wM4gAAABEhAQRAACAIUpMAQCATCtGLoqRrpLOtMUzzAwiAAAAESFBBAAAYEjNCeKaNWvitNNOi76+vsjlcnH77beXvZ8kSVx++eUxZ86c2H///WPRokXx2GOPlR3zwgsvxNKlS6OnpycOPvjgOPfcc2P79u11nQgAAMBYSrlcKh9pVHOCuGPHjjjmmGPi+uuvH/P9a6+9Nr7yla/EjTfeGOvWrYsDDzwwTjnllNi5c+fIMUuXLo2f//zncffdd8cdd9wRa9asifPPP3/fzwIAAIC61XyTmsWLF8fixYvHfC9Jkrjuuuvi05/+dLzvfe+LiIi/+7u/i1mzZsXtt98eZ555ZvziF7+Iu+66Kx588ME4/vjjIyLiq1/9aixZsiS++MUvRl9fXx2nAwAAwL5q6BrEJ554IrZu3RqLFi0aea23tzf6+/tj7dq1ERGxdu3aOPjgg0eSw4iIRYsWRT6fj3Xr1o3Z765du6JQKJQ9AAAAJqIY+VQ+0qihUW3dujUiImbNmlX2+qxZs0be27p1a8ycObPs/f322y9mzJgxckylVatWRW9v78hj7ty5jQwbAACA6JC7mK5cuTIGBgZGHlu2bGl3SAAAAJlT8xrE8cyePTsiIp599tmYM2fOyOvPPvtsHHvssSPHPPfcc2XtBgcH44UXXhhpX6m7uzu6u7sbGSoAADBJlHK5KKbsrqGZuYvpeObNmxezZ8+Oe++9d+S1QqEQ69ati4ULF0ZExMKFC2Pbtm2xfv36kWN+9KMfRalUiv7+/kaGAwAAQA1qnkHcvn17bN68eeT5E088ERs2bIgZM2bEYYcdFhdffHF89rOfjTe+8Y0xb968uOyyy6Kvry9OP/30iIg46qij4tRTT43zzjsvbrzxxtizZ09ccMEFceaZZ7qDKQAAQBvVnCA+9NBDcdJJJ408X7FiRURELFu2LG666ab41Kc+FTt27Ijzzz8/tm3bFu9617virrvuiqlTp460ufnmm+OCCy6Ik08+OfL5fJxxxhnxla98pQGnAwAAUK4U+Sil7PYraYtnWC5JkqTdQdSqUChEb29vzJ/eH125hi6jBAAAUqqYDMYjL66LgYGB6OnpqXr8cN7wv148IfbvSVfe8F+FwfiT6f9nwufSKulMWwEAAGi5dKXRAAAADVZM4V1M0xbPMDOIAAAARIQEEQAAgCFKTAEAgEwrRj6KKZsbS1s8w9IZFQAAAC0nQQQAACAilJgCAAAZp8R04tIZFQAAAC0nQQQAACAilJgCAAAZN5jrisFcV7vDKDOYS9odwpjMIAIAABAREkQAAACGKDEFAAAyrRS51N01tBS5docwpnRdJQAAANpGgggAAEBEKDEFAAAybjDyMZiyubG0xTMsnVEBAADQchJEAAAAIkKJKQAAkHHF6IpidLU7jDLFKLU7hDGZQQQAACAiJIgAAAAMUWIKABm0/ss7GtbXcRcd2LC+6uW8qnNezZfV88qyweiKwZSVmA4qMQUAACDNJIgAAABEhBJTAAAg4wYjn7qN6dMWzzAJYsrVUuPe6Br2do1da13/ZBw7K591O8fulM+60WOTXZXfq8vO2di4vlcvKHveyu+k89qHvp1Xw2X1vGAs6UxbAQAAaDkziAAAQKaVoiuKKbuLacldTAEAAEizXJIkSbuDqFWhUIje3t6YP70/unLZmgRdUljbsL7u7FnYEWM3ctzJOnanfNbtHDsrn3WtY5NdzVwTVc3VTVwz5bwaz3nVLqvnlQXFZDAeeXFdDAwMRE9PT9Xjh/OGTw38j+jueU0LIpy4XYU9cW3vbRM+l1bJVnYFAABQYTDpiq4kXSWmg4kSUwAAAFJMgggAAEBEKDFtu9E17s3ru7KGPS1jN3LcamM385wr+0/r9Z6sY0+W7xnZVv6dbt2aqEqVY9ezj1s713pVcl7VOa/ma+R58d8GIx9dKbuL6WAU2x3CmMwgAgAAEBESRAAAAIYoMQUAADJtMLpSWGKarniGmUEEAAAgIiSIAAAADFFiCgAAZFoxulJX0llMWTzDzCACAAAQERJEAAAAhigxBQAAMm0w2S/ySbpSn8Gk1O4QxmQGEQAAgIgwg9h2x110YNnzJU3sO61jN3LcamM385zbOXanfNbtHDsrn3W1sQEA6iFBBAAAMm0wuiKfsruGpu2uqsOUmAIAABAREkQAAACGKDFNmTt7FpY9X//lHRNuW7kuqZa2e9u3Z+x6xp2sY3fqZ93OsTv1sx5rbACgNoNJV+STdJV0DqYsnmFmEAEAAIgICSIAAABDlJgCAACZNhj5FN7FNJ1zdbkkSZJ2B1GrQqEQvb29MX96f3Tl5LjDlhTWNrS/yvWQrRq7lnEn69hZ+azbOXanfNbwairXsl52zsaWjX316gVlzxu5TtZ5NZ7zql1WzysLislgPPLiuhgYGIienp6qxw/nDR964eMxpae7BRFO3O7CrvjujBsnfC6tks60FQAAgJYz/QYAAGRaMdkvBpN0pT7FpNjuEMZkBhEAAICIMIPY0UbXxzev/2p7wDVy7Mq+jT36/WZ+1pVj+561fmyYiFHf2Yp1S41cM9XKNVHOq37Oq35ZPS+YCAkiAACQaYNJV+RStjH9YMriGabEFAAAgIgwgwgAmVCtJK6RfbeS86q/71ZyXvX3De0mQQQAADKtGPtFPmWpTzHcxRQAAIAUkyACAAAQEUpMASCTsrquyXl1FudFWriL6cSZQQQAACAiJIgAAAAMUWIKAABk2mCST2GJaTrn6iSIHayy/n1Jk/tv1djV6vqNnZ3Pup1jp/mzZrT1X95R0/GucWvV+vmMJ02fnfOqznlBe1xzzTWxcuXKuOiii+K6665raN/pTFsBAAAY5cEHH4yvfe1rsWDBgqb0L0EEAAAyrRj7xWDKHsV9KObcvn17LF26NP72b/82pk+f3oQrJUEEAABom0KhUPbYtWvXqx67fPnyeO973xuLFi1qWjzWIGbInT0Ly57Xu2anlvbHXbTvY9cz7mQdu55xJ+vYafqsGa3yel52zsaK57X1V7lm9OrV/12G4/OoX7XPq66+V5eXTLXy83Je+9C384K6zZ07t+z5FVdcEVdeeeWo477zne/Eww8/HA8++GBT45EgAgAAmVZMulJ3F9PiUDxbtmyJnp6ekde7u7tHHbtly5a46KKL4u67746pU6c2NS4JIgAAQJv09PSUJYhjWb9+fTz33HPxtre9beS1YrEYa9asib/+67+OXbt2RVdXYxJgCSIAAECKnXzyyfHII4+UvXbOOefEm9/85rj00ksblhxGROSSJEka1luLFAqF6O3tjfnT+6MrJ8dtlCWFtQ3rq3I9ZKvGnaxj1zLuZB27nZ811TX68xmPz652zVzrVc3VTVwL5rwaz3nRbMVkMB55cV0MDAxUnXWL+O+84YStn4n9eppbmlmrwcLO+D+zL5/wuVQ68cQT49hjj7UPIgAAAM1h+g0AAKDD3HfffU3pV4IIAABkWjHZL3JJulKfYsriGZbOqGiJ0bX8zet7vD3oGjlutbGbec6V/Vfbd69d13uyjt3K7xkT08zPZ6LjRvjsXk3559O6tV6VKseuZ3+6dq5hq+S8qsvqeUHaWYMIAABAREgQAQAAGKLEFAAAyLTBpCsiadxegY0wmLJ4hplBBAAAICKakCAeccQRkcvlRj2WL18eEXs3dKx87+Mf/3ijwwAAAKBGDS8xffDBB6NYLI48f/TRR+N3f/d34wMf+MDIa+edd1585jOfGXl+wAEHNDoMAACAiIgoJfkopqyks5Sks5iz4Qnia1/72rLn11xzTRx55JHx27/92yOvHXDAATF79uxGDw0AAEAdmpq27t69O771rW/FRz/60cjlciOv33zzzXHooYfG0UcfHStXroyXX365mWEAAAAwAU29i+ntt98e27Zti7PPPnvktQ9/+MNx+OGHR19fX2zcuDEuvfTS2LRpU9x6662v2s+uXbti165dI88LhUIzwwYAADKkmOwXkaRrA4diyuIZ1tSovv71r8fixYujr69v5LXzzz9/5N/z58+POXPmxMknnxyPP/54HHnkkWP2s2rVqrjqqquaGSoAAMCk17QE8cknn4x77rln3JnBiIj+/v6IiNi8efOrJogrV66MFStWjDwvFAoxd+7cxgU7SR130YFlz5c0se/x3m/kuNXGbuY5t3PsWq73ZB27lZ81E9PMz2ei4wIA5ZqWIK5evTpmzpwZ733ve8c9bsOGDRERMWfOnFc9pru7O7q7uxsZHgAAMEkUk66IlN3FNG13VR3WlASxVCrF6tWrY9myZbHffv89xOOPPx633HJLLFmyJA455JDYuHFjXHLJJfHud787FixY0IxQAAAAmKCmJIj33HNPPPXUU/HRj3607PUpU6bEPffcE9ddd13s2LEj5s6dG2eccUZ8+tOfbkYYAAAA1KApCeJ73vOeSJJk1Otz586N+++/vxlD0gB39iwse77+yzsm3LZyTU9tbfd93Mk6dj3jTtaxG/1Z01hXry6vIrnsnI1N6//OixraNQAdQInpxDV1H0QAAAA6hwQRAACAiGjyPogAAADtVkq6Ipeyks5SyuIZlkvGWiyYcoVCIXp7e2P+9P7oyslx02BJYW3D+qpcC2ns5o47Wceu9bOmvawZba/K69/oNaLjqVyf2sjP1nk1nvOi2YrJYDzy4roYGBiInp6eqscP5w1vevLG6OrZvwURTlyx8F+x6fCPT/hcWkWJKQAAABGhxBQAAMi4YpJP4V1M0zlXl86oAAAAaDkziOyT0esAmtd3tb3vJuPYjRy32tjNPOfK/tN0vUkXn097jfptNnHfylau9XJe9XNekD0SRAAAINNKSVfkSukqMU3rXUyVmAIAABAREkQAAACGKDEFoCE6da/CTo27naqtBWtk363kvOrvu5Wyel40RynpSt1dTJWYAgAAkGoSRAAAACJCiSkAAJBxpSQfkbKN6Uspi2eYBBGACRm9L+XGiue19bfkFf9u9p5j5fuINi7uiPLYJ+u6payet/PqLFk9L2i1dKatAAAAtJwZRAAAINOSpCuSlN01NG3xDMslSZK0O4haFQqF6O3tjfnT+6MrJ8dNgyWFtQ3r686ehcZu4biTdexaP2sa/9mPp9GfT6ti971qvlq3JRlPmkoSnVd1aTov2qeYDMYjL66LgYGB6OnpqXr8cN4wd/PNkZ92QAsinLjSSy/HljcsnfC5tIoSUwAAACJCiSkAAJBxpVJXRCldJZ2llMUzzAwiAAAAEWEGkQapXHdTy5qDyjUFta5XOO6iyTd2PeNO1rGtXand6G0t2jd2rZ9fu2KvN25Gq7a9Sl19N3l7lXHHdl61993G84LJRIIIAABkWpLkI0nZxvRpi2dYOqMCAACg5SSIAAAARIQSU5qklnUBlfuT1bpWaMmosfd9X71OGbuecSfz2ED6NXMNW6XKvpu5xs151a+V50X2lEr5iFK65sZKKYtnWDqjAgAAoOUkiAAAAESEElMAACDj9t7FNF0b06f1LqYSRNrilWsWGr0/2Xj7kDV7b7RX9l9tv79Gjl1t77XJOjaQfq1cw1ZNI9e4Oa/msyYRmiOdaSsAAAAtZwYRAADItiS/95EmaYtnSDqjAgAAoOXMIAIwpsr1O5X7ULZy7Hrbtyp2a54A6HQSRAAAINOSUj6SlG1Mn7Z4hqUzKgAAAFpOgggAAEBEKDEFYIKurthjrJH7n1X2fedFDet6VP+N3rftlX03Om4AGqTUtfeRJmmLZ4gZRAAAACJCgggAAMAQJaa0xStvBd/o28+Pd5v5Zt/6vl1jV7u1/mQdm8Yadb17FpY9Xf/lHfvcX7NLM8tib2DcEcpKATpCkk/fxvRpi2dIOqMCAACg5SSIAAAARIQSUwAAIOvcxXTCJIi03Z0NXg9US/vjLurMsesZN0tjky6d+vl0atwA0AxKTAEAAIgIM4gAAEDWJSksMU1SFs8QM4gAAABEREQuSZKk3UHUqlAoRG9vb8yf3h9dOZOgk92SwtqG9VW5HjKtYzdy3E4am/p16prPetcHt0unxt1OldfssnM2tmzsq1cvKHveyM/DeTVeM8+L9Comg/HIi+tiYGAgenp6qh4/nDfMXH9v5A9K13ektH1HPHfcyRM+l1aRXQEAANlWykWUUlY8Wcq1O4IxpewqAQAA0C4SRAAAACJCiSkdaPR6h+b1XW3Pv3aN3chxq43dzHOu7N/6kfpVWw9Uy+e3pOJ5O9cx1fq9Gy/2Rn/Pyn+bzYs7YnL8Rkb97au4Bo1c49bK6+u86jcZfw80Tq7UFbmU3cU0bfEMM4MIAABAREgQAQAAGKLEFAAAyLRc0hW5lG1Mn7Z4hkkQATKkmXuQjeq7gXtYNnvvtLL+G7z3Zqde805RbY1bI/tuJedVf99AcygxBQAAICLMIAIAAFlXyu99pEna4hmSzqgAAABoOTOIAB2s2ftU1jJ2reuDmrm350THjagv7oj2xT5Z12Nl9bydF5AWEkQAACDTcqWu1G1Mn7Z4hikxBQAAICIkiAAAAAxRYkrHqVzPsKSJfad17EaOW23sZp5ztbEB2qlyvWk90vS3LqvnBePJFbsiV0xXSWfa4hlmBhEAAICIkCACAAAwRIkpHe/OnoVlz2spnaksjam17Oa4i9ozdj3jpm1sgLQYvYXJxsb1vXpB2fNW/i3M6nlBLfKlfOTTtjF92uIZks6oAAAAaDkJIgAAABGhxBQAAMi4XCmfuo3pcyktMc0lSZK0O4haFQqF6O3tjfnT+6MrJ8dl3y0prG1of5XrIVs1di3jtntsmqvR3+nxNPKz79S4Izo79smumWvzqrm6iWv3snpeEBFRTAbjkRfXxcDAQPT09FQ9fjhvmHvvhsgfNK0FEU5caftLseXkYyd8Lq2SzrQVAACAljP9BgAAZFouyaeupDOXpCueYemMCgAAgJYzg8ikMnpdRvP6r7bXYCPHruw7TWPTWpXrfxq59qiy7zsvaljXTY27sv9Gxl3ZdzPjjmh87JNNO9fmVaocu579BLN6XkB7SBABAIBMy5e6Ip+yu5hG2uIZosQUAACAiJAgAgAAMESJKUCGjFrfU7FvXuVapVr6aub6t0bGPVZ/LYu9g+IGmExypRTexTRl8QxLZ1QAAAC0nAQRAACAiFBiCgAAZFyu2BX5YrruGpqkLJ5hEkSASaRT9yATNwC0hhJTAAAAIqIJCeKVV14ZuVyu7PHmN7955P2dO3fG8uXL45BDDomDDjoozjjjjHj22WcbHQYAAEBERORL+VQ+0qgpUb3lLW+JZ555ZuTxk5/8ZOS9Sy65JP7hH/4hvv/978f9998fTz/9dLz//e9vRhgAAADUoClrEPfbb7+YPXv2qNcHBgbi61//etxyyy3xO7/zOxERsXr16jjqqKPigQceiN/8zd9sRjgwonI90JIm99+qsautc2rn2NSunr0K26nePf/aRdwA8N+aMoP42GOPRV9fX7z+9a+PpUuXxlNPPRUREevXr489e/bEokWLRo5985vfHIcddlisXbv2VfvbtWtXFAqFsgcAAMBE5Er5VD7SqOFR9ff3x0033RR33XVX3HDDDfHEE0/Eb/3Wb8VLL70UW7dujSlTpsTBBx9c1mbWrFmxdevWV+1z1apV0dvbO/KYO3duo8MGAACY9BpeYrp48eKRfy9YsCD6+/vj8MMPj+9973ux//7771OfK1eujBUrVow8LxQKkkQAAIAGa/o+iAcffHD8+q//emzevDl+93d/N3bv3h3btm0rm0V89tlnx1yzOKy7uzu6u7ubHSqT0J09C8ue17ump7a1Y/s+dj3jNnpsald5vS87Z2PF84n3Vbme9OrVC8qeN/LzamTcEeWxNzPuiPLYGxl3RHnszYw7ojz2RsYd4bcNZFu+mI98yjamT4qTpMS00vbt2+Pxxx+POXPmxHHHHRevec1r4t577x15f9OmTfHUU0/FwoULx+kFAACAZmv4DOInP/nJOO200+Lwww+Pp59+Oq644oro6uqKs846K3p7e+Pcc8+NFStWxIwZM6KnpycuvPDCWLhwoTuYAgAAtFkuSZKkkR2eeeaZsWbNmnj++efjta99bbzrXe+Kv/iLv4gjjzwyIiJ27twZf/qnfxrf/va3Y9euXXHKKafE//pf/2vcEtNKhUIhent7Y/70/ujKNb1KFl7VksKr3323VpXlrq0at9axqV+jP7/xNPKz7dS4I1oXe6fGHTE5/w5UK5tupnaWgzeT0mWarZgMxiMvrouBgYHo6empevxw3nDU//dUdB1Y/fhWKu4oxC/+4LAJn0urNDy7+s53vjPu+1OnTo3rr78+rr/++kYPDQAAQB3SuTISAACAllOfCQAAZFq+lI98yu4ampTSFc8wCSLUYPS6jub1Pd5WFo0cdyJjU59mfm9qHbuWz7ZT4x6rfati79S4xxp7MvwdGPV3tmL9XCPX7rVybV5WzwtojXSmrQAAALScGUQAACDTcqV85FJW0pm2eIalMyoAAABazgwiAEBUX7vXyL5bKavnBTSHBBEAAMi0fDEX+WKu3WGUSVIWzzAlpgAAAESEBBEAAIAhSkwBAMaQ1fV1WT0vGE++lI98yu4amqQsnmHpjAoAAICWkyACAAAQEUpMAQCAjMuXcpEvpeuuoUnK4hkmQYQaVK7bWNLEvsd7v5HjTmRs6tPM702tY9fTtlPiHqt9q2Lv1LjHGpvWWv/lHQ3ry2cJ1EOJKQAAABFhBhEAAMi4XDEf+WK65sZKKYtnWDqjAgAAoOXMIEId7uxZWPa8ljUklWtEamu77+OONTatdfXqBWXPLztnY9P6vvOihnXdsXFX9t/IuCv7bmbcEZ11zRlf5d/tRn626ys+W3/zgVpIEAEAgEzLFXORK6brrqFpi2eYElMAAAAiIiKXJEnS7iBqVSgUore3N+ZP74+unElQOtOSwtqG9VVZ6kpnq6dUuZ06tdRZ3LRCM0tKq6ksJ/ZdoJMVk8F45MV1MTAwED09PVWPH84bjv/b52O/A6of30qDLxfiofMOmfC5tIrsCgAAyLR8kot8yjamzyfpimeYElMAAAAiQoIIAADAECWm0CKj1580r2/rSzpbp35+4m6tTo17smjnmsNKlWPbBoPJyF1MJ84MIgAAABEhQQQAAGCIElMAACDT8qUU3sU0ZfEMkyAC7IPJsFdhp8YdkZ7YOzXuiM79rgBk1apVq+LWW2+NX/7yl7H//vvHO9/5zvj85z8fb3rTmxo6jhJTAACAlLv//vtj+fLl8cADD8Tdd98de/bsife85z2xY0dt/2dkNWYQAQCATMsV9z7SpNZ47rrrrrLnN910U8ycOTPWr18f7373uxsWlwQRAACgTQqFQtnz7u7u6O7urtpuYGAgIiJmzJjR0HgkiABjqLaHWS37WC6peH51E/cgE/de48Xe6PVyr4y9mXFHpPeatzJugKyZO3du2fMrrrgirrzyynHblEqluPjii+OEE06Io48+uqHxSBABAIBMS/NdTLds2RI9PT0jr09k9nD58uXx6KOPxk9+8pOGxyVBBAAAaJOenp6yBLGaCy64IO64445Ys2ZNvO51r2t4PBJEAACAlEuSJC688MK47bbb4r777ot58+Y1ZRwJIrRI5RqcyjU7jeyb+lWux2pq3z0Lm9d3A3Vq3KP6b2Dco/pusE695s2MG2BfZOEupsuXL49bbrklfvCDH8S0adNi69atERHR29sb+++/f8Pisg8iAABAyt1www0xMDAQJ554YsyZM2fk8d3vfreh45hBBAAASLkkSVoyjgQRAADItFxp7yNN0hbPMAkitMmdFWtyKvckG481h403ek+49o1dy+fbqXFXtu/UuCPaF3unxh3hbxhAmlmDCAAAQESYQQQAADIuV0oiV2zNGr6JypXSFc8wM4gAAABEhBlESA1rcgCyo/Jv+vrVC8qeN3vvz1e6umJs/70BxiNBBAAAMs1dTCdOiSkAAAARYQYRAKDpWllyqqQUqIcEEQAAyLRcMYV3MU1ZPMOUmAIAABAREkQAAACGKDEFiNFrdJa0cex62nZK3JXtOzXuiM6JvVPjzqpqaxIb2TcQEaVk7yNN0hbPEDOIAAAARIQEEQAAgCFKTAEAgExzF9OJkyACjKFyH7Fm7lF250UN61rcE+i/kXFX9t3MuCM655o3M+6ssm4QSAslpgAAAESEGUQAACDr3MV0wswgAgAAEBERuSRJ0pm6jqNQKERvb2/Mn94fXTmToEDrrf/yjgkfm6a1RZMh7oj0xN6pcUd07ncFyLZiMhiPvLguBgYGoqenp+rxw3nD71z2n7Hf1OrHt9LgzkL86OpDJ3wurSK7AgAAsq1Y2vtIk7TFM0SJKQAAABEhQQQAAGCIElOAfdCpa67E3VqdGndEZ8cOMEqptPeRJmmLZ4gZRAAAACJCgggAAMAQJaYAAEC2lVJ4F9OUlphKEIHM6tR93CZD3BHpiV3crdep33GAyUCJKQAAABFhBhEAAMi4XKkUuZSVdKYtnmFmEAEAAIiIiFySJEm7g6hVoVCI3t7emD+9P7pyJkFhsqpcx3TZORsb1vfVqxeUPW/kOihxj9bMuCPKY29k3BHlsTcz7ojOueadGjeQfsVkMB55cV0MDAxET09P1eOH84aTP/nvsV939eNbaXBXIe794usmfC6tIrsCAACyrZjCu5imLZ4hSkwBAACICCWmQAdbUljbsrHu7FnYsL7EXV0j445oXeydGneE7wrQGfa5xPSSp9JZYvpXhykxBQAAaKlSKX0b06ctniFKTAEAAIgICSIAAABDlJgCHWP0rfPbN3Ytt9YXd+3qiXus9q2KvVPjHmvsyfJdASaJYimiWGx3FOXcxRQAAIA0kyACAAAQEUpMAQCArHMX0wkzgwgAAEBESBABAAAY0vAEcdWqVfH2t789pk2bFjNnzozTTz89Nm3aVHbMiSeeGLlcruzx8Y9/vNGhAAAA7L2DaRofKdTwBPH++++P5cuXxwMPPBB333137NmzJ97znvfEjh3lt6E+77zz4plnnhl5XHvttY0OBQAAgBo0/CY1d911V9nzm266KWbOnBnr16+Pd7/73SOvH3DAATF79uxGDw8AAMA+avoaxIGBgYiImDFjRtnrN998cxx66KFx9NFHx8qVK+Pll19+1T527doVhUKh7AEAADARSakYSTFlj1I6S0ybus1FqVSKiy++OE444YQ4+uijR17/8Ic/HIcffnj09fXFxo0b49JLL41NmzbFrbfeOmY/q1atiquuuqqZoQIAAEx6TU0Qly9fHo8++mj85Cc/KXv9/PPPH/n3/PnzY86cOXHyySfH448/HkceeeSoflauXBkrVqwYeV4oFGLu3LnNCxwAAGASalqCeMEFF8Qdd9wRa9aside97nXjHtvf3x8REZs3bx4zQezu7o7u7u6mxAl0juMuOrDs+ZI2jl1PW3HXPna97VsVe6fGPdbY9bTtlLiBSaRU3PtIk7TFM6ThCWKSJHHhhRfGbbfdFvfdd1/MmzevapsNGzZERMScOXMaHQ4AAAAT1PAEcfny5XHLLbfED37wg5g2bVps3bo1IiJ6e3tj//33j8cffzxuueWWWLJkSRxyyCGxcePGuOSSS+Ld7353LFiwoNHhAAAAMEENTxBvuOGGiIg48cQTy15fvXp1nH322TFlypS455574rrrrosdO3bE3Llz44wzzohPf/rTjQ4FAAAgnRvTpy2eIU0pMR3P3Llz4/7772/0sMAkdPXq8qqDy87Z2LS+77yoYV2LewJ9NzLuyv4bGXdl382MO6Jzrnmnxg1AC/ZBBAAAoDM0dZsLAACAtiuV0nfX0FKp3RGMKZdUqwlNoUKhEL29vTF/en905eS4wNjWf3nHhI9N063yJ0PcEemJXdyt16nfcaD9islgPPLiuhgYGIienp6qxw/nDb/zhw/HflMOakGEEze4e3v86O/eNuFzaRUlpgAAAESEElMAACDrioN7H2mStniGmEEEAAAgIswgAhnWqWuXxN1a4m69To4dIOskiAAAQKYlxWIkKduYPm3xDFNiCgAAQERIEAEAABiixBQAAMi20uDeR5qkLZ4hZhABAACICAkiAAAAQ5SYAgAA2VYaTN/G9EpMAQAASDMJIgAAABGhxBQAAMi4pFiMJJ+uks6kWGx3CGMygwgAAEBESBABAAAYosQUAADIttJg+u4amrZ4hphBBAAAICIkiAAAAAxRYgoAAGRbcTAiZXcxjWLK4hliBhEAAICIkCACAAAwRIkpAACQaUlxTyT5Pe0Oo0xSTFc8w8wgAgAAEBESRAAAAIYoMQUAADItKQ6msMTUXUwBAABIMQkiAAAAEaHEFAAAyLrS4N5HmqQtniFmEAEAAIgICSIAAABDlJgCAACZlpT2pG5j+qSUrniGmUEEAAAgIiSIAAAADFFiCgAAZFtxT0QuZSWdKSt5HWYGEQAAgIiQIAIAADBEiSkAAJBpSXEwkpSVmCbFwXaHMCYziAAAAESEBBEAAIAhSkwBAIBMS0q7Iyl1tTuMMklpd7tDGJMZRAAAACJCgggAAMAQJaYAAEC2FfdE5FKW+hTTdVfVYWYQAQAAiAgJIgAAAENSNs8KAADQWElxTyS5lN3FVIkpAAAAaSZBBAAAICKUmAIAABmXlHZHUkzX3FhS2t3uEMaUrqsEAABA20gQAQAAiAglpgAAQMYlxT2RhLuYToQZRAAAACJCgggAAMAQJaYAAECmJcXdkaRsbiwpuospAAAAKSZBBAAAICKUmAIAABmXlHZHksu1O4wySUmJKQAAACkmQQQAACAilJgCAAAZt/cupikrMXUXUwAAANJMgggAAEBEKDEFAAAyLinuSWGJ6Z52hzAmM4gAAABEhAQRAACAIUpMAQCATEtKuyKJpN1hlElK7mIKAABAikkQAQAAiAglpgAAQNYVd0WSpKvENJSYAgAAkGYSRAAAACJCiSkAAJBxSXF36kpMk9KedocwJjOIAAAARIQEEQAAgCFKTAEAgExLirsiSUrtDqOMElMAAABSTYIIAABARCgxBQAAMk6J6cSZQQQAACAiJIgAAAAMaWuCeP3118cRRxwRU6dOjf7+/vjZz37WznAAAIAMSoq795aZpuqxu92XZUxtSxC/+93vxooVK+KKK66Ihx9+OI455pg45ZRT4rnnnmtXSAAAAJNa2xLEL33pS3HeeefFOeecE7/xG78RN954YxxwwAHxjW98o10hAQAApFqzqzDbchfT3bt3x/r162PlypUjr+Xz+Vi0aFGsXbt21PG7du2KXbt2jTwvFAotiRMAAOh8pdLOyCXp2sChlAzW3Ga4CvPGG2+M/v7+uO666+KUU06JTZs2xcyZMxsSV1tmEP/zP/8zisVizJo1q+z1WbNmxdatW0cdv2rVqujt7R15zJ07t1WhAgAApEIrqjDTlUa/ipUrV8aKFStGng8MDMRhhx0WxX3IugEAgM40/L//kySpsV2xGeHUZTimyurI7u7u6O7uHnV8rVWY+6otCeKhhx4aXV1d8eyzz5a9/uyzz8bs2bNHHV95kYYv4r9uW9/cQAEAgNR56aWXore3t+pxU6ZMidmzZ8e/bn2oBVHV7qCDDhpVHXnFFVfElVdeOerY8aowf/nLXzYsprYkiFOmTInjjjsu7r333jj99NMjIqJUKsW9994bF1xwQdX2fX19sWXLlkiSJA477LDYsmVL9PT0NDlqKhUKhZg7d67r30Y+g/Zy/dvL9W8v17/9fAbt5fq3R5Ik8dJLL0VfX9+Ejp86dWo88cQTsXt3OreUSJIkcrlc2WtjzR62UttKTFesWBHLli2L448/Pt7xjnfEddddFzt27Ihzzjmnatt8Ph+ve93rRmYSe3p6/DDbyPVvP59Be7n+7eX6t5fr334+g/Zy/VtvIjOHrzR16tSYOnVqk6JpnVqrMPdV27a5+NCHPhRf/OIX4/LLL49jjz02NmzYEHfdddeoKVMAAIDJ7pVVmMOGqzAXLlzYsHHaepOaCy64YEIlpQAAAJNdPVWYE9URdzF9Nd3d3XHFFVe0vU53snL9289n0F6uf3u5/u3l+refz6C9XH/a4UMf+lD8v//3/+Lyyy+PrVu3xrHHHtvwKsxcUus9YgEAAMiktq1BBAAAIF0kiAAAAESEBBEAAIAhEkQAAAAiosMTxOuvvz6OOOKImDp1avT398fPfvazdoeUSatWrYq3v/3tMW3atJg5c2acfvrpsWnTprJjTjzxxMjlcmWPj3/8422KOFuuvPLKUdf2zW9+88j7O3fujOXLl8chhxwSBx10UJxxxhmjNlBl3x1xxBGjrn8ul4vly5dHhO9+o61ZsyZOO+206Ovri1wuF7fffnvZ+0mSxOWXXx5z5syJ/fffPxYtWhSPPfZY2TEvvPBCLF26NHp6euLggw+Oc889N7Zv397Cs+hs430Ge/bsiUsvvTTmz58fBx54YPT19cUf/uEfxtNPP13Wx1i/m2uuuabFZ9KZqv0Gzj777FHX9tRTTy07xm9g31W7/mP99yCXy8UXvvCFkWN8/+l0HZsgfve7340VK1bEFVdcEQ8//HAcc8wxccopp8Rzzz3X7tAy5/7774/ly5fHAw88EHfffXfs2bMn3vOe98SOHTvKjjvvvPPimWeeGXlce+21bYo4e97ylreUXduf/OQnI+9dcskl8Q//8A/x/e9/P+6///54+umn4/3vf38bo82WBx98sOza33333RER8YEPfGDkGN/9xtmxY0ccc8wxcf3114/5/rXXXhtf+cpX4sYbb4x169bFgQceGKecckrs3Llz5JilS5fGz3/+87j77rvjjjvuiDVr1sT555/fqlPoeON9Bi+//HI8/PDDcdlll8XDDz8ct956a2zatCl+//d/f9Sxn/nMZ8p+FxdeeGErwu941X4DERGnnnpq2bX99re/Xfa+38C+q3b9X3ndn3nmmfjGN74RuVwuzjjjjLLjfP/paEmHesc73pEsX7585HmxWEz6+vqSVatWtTGqyeG5555LIiK5//77R1777d/+7eSiiy5qX1AZdsUVVyTHHHPMmO9t27Ytec1rXpN8//vfH3ntF7/4RRIRydq1a1sU4eRy0UUXJUceeWRSKpWSJPHdb6aISG677baR56VSKZk9e3byhS98YeS1bdu2Jd3d3cm3v/3tJEmS5F//9V+TiEgefPDBkWN++MMfJrlcLvmP//iPlsWeFZWfwVh+9rOfJRGRPPnkkyOvHX744clf/dVfNTe4SWCs679s2bLkfe9736u28RtonIl8/9/3vvclv/M7v1P2mu8/na4jZxB3794d69evj0WLFo28ls/nY9GiRbF27do2RjY5DAwMRETEjBkzyl6/+eab49BDD42jjz46Vq5cGS+//HI7wsukxx57LPr6+uL1r399LF26NJ566qmIiFi/fn3s2bOn7Lfw5je/OQ477DC/hSbYvXt3fOtb34qPfvSjkcvlRl733W+NJ554IrZu3Vr2fe/t7Y3+/v6R7/vatWvj4IMPjuOPP37kmEWLFkU+n49169a1PObJYGBgIHK5XBx88MFlr19zzTVxyCGHxFvf+tb4whe+EIODg+0JMIPuu+++mDlzZrzpTW+KP/7jP47nn39+5D2/gdZ59tln4x//8R/j3HPPHfWe7z+dbL92B7Av/vM//zOKxWLMmjWr7PVZs2bFL3/5yzZFNTmUSqW4+OKL44QTToijjz565PUPf/jDcfjhh0dfX19s3LgxLr300ti0aVPceuutbYw2G/r7++Omm26KN73pTfHMM8/EVVddFb/1W78Vjz76aGzdujWmTJky6n+YzZo1K7Zu3dqegDPs9ttvj23btsXZZ5898prvfusMf6fH+ts//N7WrVtj5syZZe/vt99+MWPGDL+JJti5c2dceumlcdZZZ0VPT8/I65/4xCfibW97W8yYMSN++tOfxsqVK+OZZ56JL33pS22MNhtOPfXUeP/73x/z5s2Lxx9/PP78z/88Fi9eHGvXro2uri6/gRb65je/GdOmTRu1rMP3n07XkQki7bN8+fJ49NFHy9bARUTZ2ob58+fHnDlz4uSTT47HH388jjzyyFaHmSmLFy8e+feCBQuiv78/Dj/88Pje974X+++/fxsjm3y+/vWvx+LFi6Ovr2/kNd99Jqs9e/bEBz/4wUiSJG644Yay91asWDHy7wULFsSUKVPiYx/7WKxatSq6u7tbHWqmnHnmmSP/nj9/fixYsCCOPPLIuO++++Lkk09uY2STzze+8Y1YunRpTJ06tex13386XUeWmB566KHR1dU16k6Nzz77bMyePbtNUWXfBRdcEHfccUf8+Mc/jte97nXjHtvf3x8REZs3b25FaJPKwQcfHL/+678emzdvjtmzZ8fu3btj27ZtZcf4LTTek08+Gffcc0/80R/90bjH+e43z/B3ery//bNnzx51s7LBwcF44YUX/CYaaDg5fPLJJ+Puu+8umz0cS39/fwwODsavfvWr1gQ4ibz+9a+PQw89dORvjt9Aa/zzP/9zbNq0qep/EyJ8/+k8HZkgTpkyJY477ri49957R14rlUpx7733xsKFC9sYWTYlSRIXXHBB3HbbbfGjH/0o5s2bV7XNhg0bIiJizpw5TY5u8tm+fXs8/vjjMWfOnDjuuOPiNa95TdlvYdOmTfHUU0/5LTTY6tWrY+bMmfHe97533ON895tn3rx5MXv27LLve6FQiHXr1o183xcuXBjbtm2L9evXjxzzox/9KEql0kjyTn2Gk8PHHnss7rnnnjjkkEOqttmwYUPk8/lRpY/U79///d/j+eefH/mb4zfQGl//+tfjuOOOi2OOOabqsb7/dJqOLTFdsWJFLFu2LI4//vh4xzveEdddd13s2LEjzjnnnHaHljnLly+PW265JX7wgx/EtGnTRtYw9Pb2xv777x+PP/543HLLLbFkyZI45JBDYuPGjXHJJZfEu9/97liwYEGbo+98n/zkJ+O0006Lww8/PJ5++um44ooroqurK84666zo7e2Nc889N1asWBEzZsyInp6euPDCC2PhwoXxm7/5m+0OPTNKpVKsXr06li1bFvvt999/Nn33G2/79u1ls69PPPFEbNiwIWbMmBGHHXZYXHzxxfHZz3423vjGN8a8efPisssui76+vjj99NMjIuKoo46KU089Nc4777y48cYbY8+ePXHBBRfEmWeeWVYazKsb7zOYM2dO/MEf/EE8/PDDcccdd0SxWBz5b8KMGTNiypQpsXbt2li3bl2cdNJJMW3atFi7dm1ccskl8ZGPfCSmT5/ertPqGONd/xkzZsRVV10VZ5xxRsyePTsef/zx+NSnPhVveMMb4pRTTokIv4F6VfsbFLH3/5j6/ve/H3/5l385qr3vP5nQ7tuo1uOrX/1qcthhhyVTpkxJ3vGOdyQPPPBAu0PKpIgY87F69eokSZLkqaeeSt797ncnM2bMSLq7u5M3vOENyZ/92Z8lAwMD7Q08Iz70oQ8lc+bMSaZMmZL82q/9WvKhD30o2bx588j7//Vf/5X8yZ/8STJ9+vTkgAMOSP7H//gfyTPPPNPGiLPnn/7pn5KISDZt2lT2uu9+4/34xz8e8+/NsmXLkiTZu9XFZZddlsyaNSvp7u5OTj755FGfy/PPP5+cddZZyUEHHZT09PQk55xzTvLSSy+14Ww603ifwRNPPPGq/0348Y9/nCRJkqxfvz7p7+9Pent7k6lTpyZHHXVU8rnPfS7ZuXNne0+sQ4x3/V9++eXkPe95T/La1742ec1rXpMcfvjhyXnnnZds3bq1rA+/gX1X7W9QkiTJ1772tWT//fdPtm3bNqq97z9ZkEuSJGl6FgoAAEDqdeQaRAAAABpPgggAAEBESBABAAAYIkEEAAAgIiSIAAAADJEgAgAAEBESRAAAAIZIEAEAAIgICSIAAABDJIgAAABEhAQRAACAIRJEAAAAIiLi/w+nVuRpY6UscAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x1000 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(16, 10), dpi=100)\n",
    "norm = plt.Normalize(0, 10)\n",
    "cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\", \"orange\"])\n",
    "cbar = fig.colorbar(\n",
    "    ax.imshow(phantom.image.T, origin=\"lower\", cmap=mpl.colormaps['turbo'], norm=norm, aspect=\"equal\"), pad=0.01\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = np.zeros((180, 180))\n",
    "mask = np.zeros((180, 180))\n",
    "radii = [2,4,6,8,10,12]\n",
    "R=55\n",
    "for idx in range(0,6):\n",
    "    angle_rad = idx*math.pi/3\n",
    "    cx,cy = (math.cos(angle_rad)*R+90,math.sin(angle_rad)*R+90)\n",
    "    xx,yy = skimage.draw.circle_perimeter(int(cx),int(cy),radii[idx])\n",
    "    img[xx,yy] = 5\n",
    "    mask[xx,yy] = idx+1\n",
    "    xx,yy = skimage.draw.disk((int(cx),int(cy)),radii[idx])\n",
    "    img[xx,yy] = 10\n",
    "    mask[xx,yy] = idx+1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(16,10),dpi=100)\n",
    "norm=plt.Normalize(0,10)\n",
    "cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\"])\n",
    "fig.colorbar(ax.imshow(img.T,origin='lower',cmap=cmap,norm=norm))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataDir = \"/home/fanghan/Work/SPEBT/Develop/data/2023/\"\n",
    "NImgX_ = 180\n",
    "NImgY_ = 180\n",
    "NDetY_ = 32\n",
    "NModule_ = 6\n",
    "NDetX_ = 24\n",
    "fName = (\n",
    "    \"sysmat_24layer_Rot_0_of_2_2mmslitin10mm_1010_idxT0_numT1in1mm_IZ0_DZ0_100.sysmat\"\n",
    ")\n",
    "data = np.fromfile(\"%s%s\" % (dataDir, fName), dtype=\"single\")\n",
    "dataArr2D= data.reshape((NDetX_ * NModule_ * NDetY_, NImgX_ * NImgY_))\n",
    "dataArr5D = dataArr2D.reshape((NDetX_, NModule_, NDetY_, NImgX_, NImgY_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sys.getsizeof(data)/1024/1024,\"MB\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nlayers = 1\n",
    "sysmat = np.reshape(dataArr5D[:Nlayers,:,:,:,:],(Nlayers,NModule_,NDetY_,NImgX_ * NImgY_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 500 Bq/uL, 30s\n",
    "# 10 * 1500 Count/uL\n",
    "\n",
    "poisson_img = np.random.poisson(img*1500)\n",
    "arr_dict={}\n",
    "for Nlayers in range(1, 11):\n",
    "    dataMat = np.reshape(dataArr5D[:Nlayers, :, :, :, :],(Nlayers*NModule_*NDetY_,NImgX_ * NImgY_))\n",
    "    sysmat = dataMat[~np.all(dataMat == 0, axis=1)]\n",
    "    fake_proj = np.matmul(sysmat, poisson_img.reshape(NImgX_ * NImgY_))\n",
    "    arr_dict['%d-layers'%Nlayers] = fake_proj\n",
    "    print(fake_proj.shape)\n",
    "# for keys in arr_dict:\n",
    "#     print(arr_dict[keys].shape)\n",
    "np.savez(\"test.npz\",**arr_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npzfile = np.load(\"test.npz\")\n",
    "sorted(npzfile.files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nlayers=8\n",
    "dataMat = np.reshape(dataArr5D[:Nlayers, :, :, :, :],(Nlayers*NModule_*NDetY_,NImgX_ * NImgY_))\n",
    "# print(dataMat.shape)\n",
    "sysmat = dataMat[~np.all(dataMat == 0, axis=1)]\n",
    "print(sysmat.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def backwardProj(lastArr, projArr, sysMat):\n",
    "    forwardLast = np.matmul(sysMat, lastArr)\n",
    "    quotients = projArr/forwardLast\n",
    "    return np.matmul(quotients, sysMat)/np.sum(sysMat, axis=0)*lastArr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "Nlayers=8\n",
    "dataMatrix = np.reshape(dataArr5D[:Nlayers, :, :, :, :],(Nlayers*NModule_*NDetY_,NImgX_ * NImgY_))\n",
    "sysMatrix = dataMatrix[~np.all(dataMatrix == 0, axis=1)]\n",
    "matSum = np.sum(sysMatrix, axis=0, dtype=float)\n",
    "validIdx = np.array(np.where(matSum != 0))\n",
    "# validIdx=np.where(matSum != 0)\n",
    "reducedMat = np.reshape(\n",
    "    sysMatrix[:, validIdx], (sysMatrix.shape[0], validIdx.shape[1]))\n",
    "projection = np.random.poisson(arr_dict['%d-layers'%Nlayers])\n",
    "print(projection.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = projection == np.floor(projection)\n",
    "np.alltrue(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poisson_img = np.random.poisson(img*1500)\n",
    "fig,ax = plt.subplots(figsize=(16,10),dpi=100)\n",
    "# norm=plt.Normalize(0,10)\n",
    "# cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\"])\n",
    "fig.colorbar(ax.imshow(poisson_img.T,origin='lower'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def myfunc(k,mean):\n",
    "    return mean**k * math.exp(-mean)/math.factorial(k)\n",
    "vmyfunc = np.vectorize(myfunc)\n",
    "\n",
    "def my_log_gaus(xx,mean):\n",
    "    return np.log(1/np.sqrt(2*3.1416*mean)*np.exp(-(xx-mean)**2/(2*mean)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(projection),np.std(projection))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NIteration = 50\n",
    "print((reducedMat.shape))\n",
    "reconImg = np.ones(reducedMat.shape[1])*200\n",
    "for iter in range(NIteration):\n",
    "        reconImg = backwardProj(reconImg, projection, reducedMat)\n",
    "        probs = my_log_gaus(projection,np.matmul(reducedMat,reconImg))\n",
    "        print(np.sum(probs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = np.arange(1,5)\n",
    "means = np.array([0.5,2.1,3.6,4.7])\n",
    "vmyfunc = np.vectorize(myfunc)\n",
    "print(vmyfunc(ks,means))\n",
    "for k,mean in zip(ks,means):\n",
    "    print(myfunc(k,mean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(16,10),dpi=100)\n",
    "plotImg = np.zeros((NImgX_ * NImgY_))\n",
    "plotImg[validIdx] = reconImg\n",
    "# norm=plt.Normalize(0,10)\n",
    "# cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\"])\n",
    "fig.colorbar(ax.imshow(np.reshape(plotImg,(NImgX_,NImgY_)).T,origin='lower'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "NIteration = 100\n",
    "Nlayers=2\n",
    "dataMatrix = np.reshape(dataArr5D[:Nlayers, :, :, :, :],(Nlayers*NModule_*NDetY_,NImgX_ * NImgY_))\n",
    "sysMatrix = dataMatrix[~np.all(dataMatrix == 0, axis=1)]\n",
    "matSum = np.sum(sysMatrix, axis=0, dtype=float)\n",
    "validIdx = np.array(np.where(matSum != 0))\n",
    "# validIdx=np.where(matSum != 0)\n",
    "reducedMat = np.reshape(\n",
    "    sysMatrix[:, validIdx], (sysMatrix.shape[0], validIdx.shape[1]))\n",
    "projection = np.random.poisson(arr_dict['%d-layers'%Nlayers])\n",
    "\n",
    "reconImg = np.ones(reducedMat.shape[1])*200\n",
    "for iter in range(NIteration):\n",
    "        reconImg = backwardProj(reconImg, projection, reducedMat)\n",
    "        # probs = my_log_gaus(projection,np.matmul(reducedMat,reconImg))\n",
    "        # print(np.sum(probs))\n",
    "fig,ax = plt.subplots(figsize=(16,10),dpi=100)\n",
    "plotImg = np.zeros((NImgX_ * NImgY_))\n",
    "plotImg[validIdx] = reconImg\n",
    "# norm=plt.Normalize(0,10)\n",
    "# cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\",\"orange\"])\n",
    "fig.colorbar(ax.imshow(np.reshape(plotImg,(NImgX_,NImgY_)).T,origin='lower'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NIteration = 100\n",
    "for Nlayers in range(1,10):\n",
    "    dataMatrix = np.reshape(dataArr5D[:Nlayers, :, :, :, :],(Nlayers*NModule_*NDetY_,NImgX_ * NImgY_))\n",
    "    sysMatrix = dataMatrix[~np.all(dataMatrix == 0, axis=1)]\n",
    "    matSum = np.sum(sysMatrix, axis=0, dtype=float)\n",
    "    validIdx = np.array(np.where(matSum != 0))\n",
    "    reducedMat = np.reshape(\n",
    "        sysMatrix[:, validIdx], (sysMatrix.shape[0], validIdx.shape[1]))\n",
    "    projection = np.random.poisson(arr_dict['%d-layers'%Nlayers])\n",
    "\n",
    "    reconImg = np.ones(reducedMat.shape[1])*200\n",
    "    for iter in range(NIteration):\n",
    "            reconImg = backwardProj(reconImg, projection, reducedMat)\n",
    "    fig,ax = plt.subplots(figsize=(16,10),dpi=100)\n",
    "    plotImg = np.zeros((NImgX_ * NImgY_))\n",
    "    plotImg[validIdx] = reconImg\n",
    "    mytext=\"N x-layers: %2d\"% Nlayers\n",
    "    ax.annotate(mytext,\n",
    "                xy=(2, 175), xycoords='data',fontsize=18,\n",
    "                horizontalalignment='left',\n",
    "                verticalalignment='center',color=\"white\")\n",
    "    ax.tick_params(labelsize=12)\n",
    "    fig.colorbar(ax.imshow(np.reshape(plotImg,(NImgX_,NImgY_)).T,origin='lower'),pad=0.01)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(\"%d-x-layers-recon.png\"%Nlayers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = np.zeros((imgN_x_,imgN_y_))\n",
    "img[89,89]=10\n",
    "fig, ax = plt.subplots(figsize=(11, 10), dpi=100)\n",
    "norm = plt.Normalize(0, 10)\n",
    "cmap = mpl.colors.LinearSegmentedColormap.from_list(\"\", [\"white\", \"orange\"])\n",
    "cbar = fig.colorbar(\n",
    "    ax.imshow(\n",
    "        img.T, origin=\"lower\", cmap=mpl.colormaps[\"gray\"], norm=norm, aspect=\"equal\"\n",
    "    ),\n",
    "    pad=0.01,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "outDir='output/'\n",
    "if os.path.exists(outDir)!= True:\n",
    "    try:\n",
    "        os.mkdir(outDir)\n",
    "    except Exception as err:\n",
    "        print(err)\n",
    "        exit(2)\n",
    "outImgName='dot_phantom_plot_%dx%d.png'%(imgN_x_,imgN_y_)\n",
    "outNpzName='dot_phantom_data_%dx%d.npz'%(imgN_x_,imgN_y_)\n",
    "print(\"Save phantom plot to\",outDir+outImgName)\n",
    "fig.savefig(outDir+outImgName,dpi=100)\n",
    "print(\"Save phantom data to\",outDir+outNpzName)\n",
    "np.savez_compressed(outDir+outNpzName,phantom=img,mask=mask)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}