Skip to content

Instantly share code, notes, and snippets.

@pangyuteng
Last active November 16, 2022 18:34
Show Gist options
  • Save pangyuteng/d062252d4eb02e25d9ea4b7ffd5bcd54 to your computer and use it in GitHub Desktop.
Save pangyuteng/d062252d4eb02e25d9ea4b7ffd5bcd54 to your computer and use it in GitHub Desktop.
3d image patch extraction
**/.ipynb_checkpoints/*
ok.png
ok.nii.gz
workdir/*

problem statement

When moving the patches generated by tf.extract_volume_patches or tf.image.extract_patches back to its original shape with tf.reshape or np.reshape methods, the result will be incorrect.

Above issue is explained in below extensively.

https://stackoverflow.com/questions/44047753/reconstructing-an-image-after-using-extract-image-patches

Note Stop right here if just need non overlapping patches:

  • see if np.reshape fits your bill - which usually does, then use that! (np-extract-patches.ipynb)

https://stackoverflow.com/questions/31527755/extract-blocks-or-patches-from-numpy-array

  • also for 2d images, for non-overlapping patches, one can use tf.nn.space_to_depth and tf.nn.depth_to_space. See tf-space-to-depth.ipynb.

  • alternatively if you prefer c++, git ITK a try.

https://discourse.itk.org/t/patch-extraction/3630

https://itk.org/ITKExamples/src/Filtering/ImageGrid/ExtractRegionOfInterestInOneImage/Documentation.html

solutions

Several solutions to resolve the above issue are listed below:

notebooks can be executed following below prep script within tf container.
docker run -it -p 8888:8888 -w /workdir -v $PWD:/workdir tensorflow/tensorflow:2.6.1-gpu-jupyter bash
bash prep.sh
jupyter notebook --ip=* --allow-root

additional thoughts

whatever solutions you end up using, doesn't really matter during training.

You will for sure get a performance hit during inference. For problems dealing with 3d patches, tf.extract_volume_patches may be one method, but for medical images, depending on hardware, OOM may occur.

  • thus now you are back to for loops with python/numpy for preprocessing -> batch wise inference -> reshape

  • alternatively, chop up the volume to cubes (faster than for loops) -> feed that asyncronously to multiple serving end points -> gather the results and do the reshape.

  • wait until you get the budget to run jumbo gpus for inference 24/7.

reference


https://stackoverflow.com/questions/41564321/split-image-tensor-into-small-patches
https://stackoverflow.com/questions/44047753/reconstructing-an-image-after-using-extract-image-patches
https://stackoverflow.com/questions/40731433/understanding-tf-extract-image-patches-for-extracting-patches-from-an-image

https://stackoverflow.com/questions/53451728/n-dimensional-sliding-window-operation-in-python-using-gpu-accelerated-libraries
https://stackoverflow.com/questions/35691947/3d-sliding-window-operation-in-theano


vae patch
Discovering Digital Tumor Signatures—Using Latent Code Representations to Manipulate and Classify Liver Lesions
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8269051

pediatric-ct-seg
https://wiki.cancerimagingarchive.net/pages/viewpage.action?pageId=89096588#89096588bcab02c187174a288dbcbf95d26179e8

Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import numpy as np\n",
"import SimpleITK as sitk\n",
"from scipy.ndimage import binary_dilation"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"\n",
"def imread(fpath):\n",
" if fpath.endswith('.list'):\n",
" with open(fpath,'r') as f:\n",
" dicom_names = [x for x in f.read().split('\\n') if len(x) > 0]\n",
" if not os.path.exists(dicom_names[0]):\n",
" dicom_names = [os.path.join(os.path.dirname(fpath),x) for x in dicom_names]\n",
" reader = sitk.ImageSeriesReader()\n",
" reader.SetFileNames(dicom_names)\n",
" else:\n",
" reader= sitk.ImageFileReader()\n",
" reader.SetFileName(fpath)\n",
" img = reader.Execute()\n",
" #arr = sitk.GetArrayFromImage(img)\n",
" #spacing = img.GetSpacing()\n",
" #origin = img.GetOrigin()\n",
" #direction = img.GetDirection() \n",
" return img\n",
"\n",
"def imwrite(fpath,arr,spacing,origin,direction,use_compression=True):\n",
" img = sitk.GetImageFromArray(arr)\n",
" img.SetSpacing(spacing)\n",
" img.SetOrigin(origin)\n",
" img.SetDirection(direction)\n",
" writer = sitk.ImageFileWriter() \n",
" writer.SetFileName(fpath)\n",
" writer.SetUseCompression(use_compression)\n",
" writer.Execute(img)\n",
" \n",
"\n",
"# ref https://gist.github.com/mrajchl/ccbd5ed12eb68e0c1afc5da116af614a\n",
"def resample_img(itk_image, out_spacing=[1.0, 1.0, 1.0], is_label=False):\n",
" \n",
" # Resample images to out_spacing with SimpleITK\n",
" original_spacing = itk_image.GetSpacing()\n",
" original_size = itk_image.GetSize()\n",
"\n",
" out_size = [\n",
" int(np.round(original_size[0] * (original_spacing[0] / out_spacing[0]))),\n",
" int(np.round(original_size[1] * (original_spacing[1] / out_spacing[1]))),\n",
" int(np.round(original_size[2] * (original_spacing[2] / out_spacing[2])))]\n",
"\n",
" resample = sitk.ResampleImageFilter()\n",
" resample.SetOutputSpacing(out_spacing)\n",
" resample.SetSize(out_size)\n",
" resample.SetOutputDirection(itk_image.GetDirection())\n",
" resample.SetOutputOrigin(itk_image.GetOrigin())\n",
" resample.SetTransform(sitk.Transform())\n",
" resample.SetDefaultPixelValue(itk_image.GetPixelIDValue())\n",
"\n",
" if is_label:\n",
" resample.SetInterpolator(sitk.sitkNearestNeighbor)\n",
" else:\n",
" resample.SetInterpolator(sitk.sitkBSpline)\n",
"\n",
" return resample.Execute(itk_image)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['rlung.nii.gz',\n",
" 'wbone.nii.gz',\n",
" 'rll.nii.gz',\n",
" 'scanner.nii.gz',\n",
" 'rul.nii.gz',\n",
" 'rml.nii.gz',\n",
" 'lul.nii.gz',\n",
" 'blanket-bed.nii.gz',\n",
" 'lll.nii.gz',\n",
" 'image.nii.gz',\n",
" 'bottom-air.nii.gz',\n",
" 'llung.nii.gz',\n",
" 'top-air.nii.gz']"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"os.listdir('workdir')"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {},
"outputs": [],
"source": [
"img = resample_img(imread('workdir/image.nii.gz'))\n",
"lul = resample_img(imread('workdir/lul.nii.gz'),is_label=True)\n",
"lll = resample_img(imread('workdir/lll.nii.gz'),is_label=True)\n",
"llung = resample_img(imread('workdir/llung.nii.gz'),is_label=True)"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [],
"source": [
"img = sitk.GetArrayFromImage(img).copy()\n",
"lul = sitk.GetArrayFromImage(lul).copy()\n",
"lll = sitk.GetArrayFromImage(lll).copy()\n",
"llung = sitk.GetArrayFromImage(llung).copy()"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"((301, 312, 312), (301, 312, 312), (301, 312, 312), (301, 312, 312))"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"img.shape,lul.shape,lll.shape,llung.shape"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [],
"source": [
"# assume complete lof \n",
"overlap = binary_dilation(lul).astype(int)+binary_dilation(lll).astype(int)\n",
"fiss = np.zeros(overlap.shape)\n",
"fiss[ovelap==2]=1\n",
"fiss[llung==0]=0\n",
"fiss = (fiss > 0).astype(int)\n"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7feac70f1b50>"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAACCCAYAAAC96IjgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmlUlEQVR4nO3dW4xcyX3f8e+/qs6lu6fnxuGdXGqlXVuWb7Ks2DLi3CwLkfUiPxiO/WALhoB9cQAbyEOEBIhf8mDnwQYMBAYE2IgUGDEE2YGVRIEtKYrlJJIgyRdZ2pW0u9Zyl1ySS3I4t76cc6rqn4dzZjjcXS0pci6cZn2AQfecbnLO4Z/9m+rquoiqkiRJkswWc9gnkCRJkuy9FO5JkiQzKIV7kiTJDErhniRJMoNSuCdJksygFO5JkiQzaF/CXUTeKyLfFJHnRORD+/EzksORajubUl1nj+z1OHcRscC3gPcAl4AvAb+oqk/v6Q9KDlyq7WxKdZ1N+9Fy/zHgOVX9e1WtgT8C3r8PPyc5eKm2synVdQbtR7ifBV7a9f2l7lhy9KXazqZU1xnkDusHi8hTwFMAFvujfeYP61SSXTa5dUNVj9/vn091fTilus6mN6rrfoT7ZeD8ru/PdcfuoKofBj4MMC/L+uPy7n04leS79Wn9+MU3ePiutU11fTilus6mN6rrfnTLfAl4UkQeF5Ec+AXgE/vwc5KDl2o7m1JdZ9Cet9xV1YvIvwT+DLDAH6jq1/f65yQHL9V2NqW6zqZ96XNX1U8Cn9yPvzs5XKm2synVdfakGapJkiQzKIV7kiTJDErhniRJMoNSuCdJksygFO5JkiQz6NBmqD7SjMX0SqTfB0A3N4nT6SGfVJIksySF+wEzwyE8fha/1CNagwkRO/bYtRF69Tpxc/OwTzFJkhmQwv0AiXPI+dNUKwP8wBKdIKrIMMPN52SLA9yVVcL1G2hVHfbpJklyhKVwP0D29CmqE3M0Q4svDdEBIqBQDw12yZEdKyluLGEvXiNcv37Yp5wkyRGVwv2A2KUl/Okl6kVHPWcImaAW1IAoEMFHoekZfH9ALz9DVhbEm6vE8fiwTz9JkiMmhfsBMGWJnj1Jdaykmrc0fYi5oAYQULkd8CaA7xl80aPfO0nxcg975RXC+gbs8a5ZSZLMrhTu+81Y5MI5pqcGTI9Z6qHg+xCzttW+3XoHkACmAeOFUAq+l1MsLVOeGlJ8+wb+4ksp4JMkuScp3PeZGfSJ8z2aocOXQswh5BBKRTOITkHa54qXNtzrrlWvghpDyHN8eZIB4F+8DDEc5iUlSXIEpHDfZ2ZxgWouIxRCyAVfQugpvq9orqjb1RKPII3BVoCRtqtGBASidegPnmLOB/yl1+x9kiRJcocU7vtInCMuzOF7Fl8KoXc72GM/gotgQIyCKBoFzSO+MERniE5wTrpbRY3FvvUUZVWnkTRJkryhFO77yPT7xF5GKA2+B75PG+xlBKtgFeMixioiSlRpA94rEfBiUAfRCWrbrprJiYzs8VPYqiJsbBz2JSZJ8pBK4b6feiWhaFvtzUBo5pRYKmQKLmJcxLqIywLGtN0zMQreG4KNRGeJpSFmhpi1LXiJwvhsn+HoNPJclSY7JUnyulK47yNxjlhYmp7g5yAWXR9712K3LlKUDbnzZDa2LXcVam9pXMS7iK8cwShqbduFEwVTG9y5If3JacKLl1DvD/tSkyR5yKRw308ihNx0H6Z2we4iJg+4LFAUDfNlxSCrEWlb7qpCHS1T75jWGRMbaaxru2nUYhqhmRemxxxuukRe1fjLLx/udSZJ8tBJ4b6PtKpRK4SyG8+eRWwv0OtX9PKGpXLCQjGhtA2FuT280ath5HM26pKtomBjWjA2JUEFXwmmFqoFAS2YCydwG5tpwbEkSe6Qwn0fxVu36F07xejUHJMTYPqe+eGY44MRw3zKQjZlmE0pjCeTQCaBiNBEy5YrGLiam3aw8/eNFHwEE9oWvEQYnykZTh5D/vabqXsmSZIdKdz3kXqPe+kVlvoOtSWrxxxLpyecG6yxmI0pjGfBTShMg0XJxBMwVDHbCfygQtR2llMIhok3NF6QsL12AdjzA+ZeXErDI5Mk2ZHCfZ/FjU3sRs3c5Yx6IWPjdMnplXWW3YiAsGxHZOIZmJpMPFPNGMcCIxGAKmZENfhoqLyjaSxNFJrYLkojQZjWlv65E3BzNc1eTZIEuIdt9kTkD0TkFRH52q5jyyLyKRF5trtd6o6LiPyuiDwnIl8VkXfs58kfCSJIjNg6km0oN1fniCos2DHLdsRxt8HZ7Bbns5ucsJssmjFDM2XRjunbioGr6NmG0noK5ykKjy09sYyEUtvx8z2oTvSwc4O7n0/n6/pl/kL/G5/XP9851mgN8GSq694w/T52cQGMPdCf+3q1BWx6ze4RYxH38LeL72UP1f8EvPdVxz4EfEZVnwQ+030P8DPAk93XU8Dv7c1pHnGqmCbiJsBGxvV6yKIdc9xttLd2xEA8C6Zi3kwZmgl9UzE0U/qmpmdrStfQyxp6eUOee6TnCYNI6Cu+JzRDiyzM3/MpneECP8JP3nHsBb4BsJnqugdEkOEcMhxiBv0D/dGvV1vgNOk1uydMniFFcdincVd3DXdV/Ryw+qrD7wc+0t3/CPCzu45/VFtfABZF5PQeneuRJUGRqLipQtvbwtBMOOvWGEjNQDwGxYoyNDWD7ZC3UxbcmDlbMXQVw2zKIK8p8waXB7QMhFIJJdQDQ1hZaNeiuQdLcpyM/I5j13kZ4Gb3barrA5A8R7IMYsQszN+uyz3W50G8Xm2BRdJrds+Iffhb7/d7didV9Up3/ypwsrt/Fnhp1/Mudceu8CjTNtTVgA4C3zO4yim3yaLxrIaMTCCTSFBAIqV4gmlnnjZqqVxGFR1eDXV0TL2jyj2+scRSb3fPLBa4PL/vWas1FUDTfZvq+gC0qkAELXPEB9zjF2A8aec+rN46jJnFLr1m94Z6D8a04S6328fa1Id4Vq/1wL96VFVlewbOd0FEnqJ9G0jJwb5tPUhiLZrZnSUEzp5d5T2DpzllAxmGxjQYwAKZQOzWa7dE+qbiuNsgqMFIxJlIEy1NsPhg2w9XnSV2E6WqxYx8OEfYg+BIdX1AIpB1L6+qhjwDa8EY7NIi/uq1Qzu1+6ltqusuYpA8gyxvV+vWCHVDCOGhGtBwv+F+TUROq+qV7i3cK93xy8D5Xc871x17DVX9MPBhgHlZnt0dKEy7JjvSttxP9DfpG08hhkYjy8awGttbgLXYttYBcgLHzIjg2seqmLGYTboWvGXaOJo8IzoldOvXsLzYjpq5j009cgo8TQaQ6vpgTFGAKjKagnOE5XliP8NUHjOqkBs3D3pegn+Q12yqa0ekDXbnEGfRuT5kDhlNMN4/VFti3ssHqq/nE8AHuvsfAP501/Ff7j6Bfxewvuut4KMpKhIjdC+H51dX+KvpOS564aLPeDkIU7WYri92PRZMNQPASCST0PW/T1iwExazMUPXLlmQ2YDJIjGDmNMuK7w8QFx2X6d6nDMAx7pvU10fgCzMo3mGLs0TTiwQ5nLctXXszU38Uv8wPpBbI71mH5hYi2SufRfmHHG+h1/ooUX+0H3IeteWu4j8F+CfAisicgn4DeA3gY+JyAeBi8DPd0//JPA+4DlgDPzKPpzzkaKqiI9IbNN9c63Pn936Ab7ZO82CG9M3NW8vLzKOFQ1toPelIieyFnsMzZShmQLQqGv74KNj0xf084YtF/BOiZkScqFeLOjNDQi33rj/7+/0i9ziOg0Vf6n/gzfzNi7wvVzkW/Mi8iyprvdNigKOLRL7OeIjoZchISI+oP0Sd2sMp0/A8y/sy7aJr1db2j7096TX7IORomh/cZcFqopag4QI1iCDPmYyIU6nh32awD2Eu6r+4nd46N2v81wFfvVBT2qW6LRCGt8OhZwq5kbG8+sr9GxDJoG+qXmhWWEcN+mbikYtAcP1WHLMjBh2s1cDwsBU9E1FYTy58WQ24FyksYqadl/WUApSvGakxGv8oPz4dzhhvqWq77zjUKrrPTP9Prz5MWJm8XM5aoRQGPovboAxxCLDjKZgLabX25e38a9X22f0K0FV02v2AYhzSL+P9kvCsCQWjlAYbBORkGMnFdLrQVU9FHsdP9xjeWaBRlBFBdw0Utx0vHxzgZXeFmt1jyfnXuHZcJIT+Qbns3bEqZFIKQ1DU2NRCoFSAqU0O2vQFCbQcw1F5plk3drvOYTcwEM+RGuWiHOYpSVYHBIX+ozO9hmfsExWBDeF+RcC/SsTwrDAn5zD1AFT+3aI5PLSQ9VHm7yKSPvhqZF2aOubzhEHBdWxkmZgdja2d1MlM4IZ14gPD03rPaXAPtMQkLrB1IF8zTP/gmHd9fm6Pc3KQhvwt6Y9Suf54eXLnM7XySRwIb+B6VrslSpD8QzNhKEtGdopA1e1q0lmHpMHNHPETPBlO/wuORhmcQE9eYwwXxB6jtFJS7Us1EvK4GmQqPhBRvHyBnZ9wuTCIlghf/4VdH6AGQyIo9FhX0byKuIc9uQJdHFILB2+lzE+XTBdNtTDdl9jW4GtFNNAUQhus8D40E4mfAha7ync95squr6Bnetj6sBcUHxZcnOlzzVvWe/1GK2XiFFCNHzfUjsiJrSDrDjl1sgkcMxULJqKta5rpm9qSutxpt30w7uuW6YA7qFbJtkjJ45Rnxzg+xaJyviMkK/Bmc8FetcmqBXqhZzRW5YwTeTm92eMzllOfOk82SiSn5jDff7raUeth4xZmCecWCLM5fi+ZbLi2Lxg2o3tbftO3DTgxoKtwPcsauaYe9FgjcFOpoe+DWYK9wMQ1taRrRGIwRrh+JUz+MFpNrRg1HcQ2o2x1yYlz3CK5XLMRbfMN+xpHu/fYMVtcia7xXz3wWppGvq2bvvtbcDYiOaR6ISQCZod7FomjzKZVJh6gOSGmAun/19N8fLW7ceBIiihdKgTTv/fLcZnelz/YUvvumHrJyzuH/8op79QU357FW6tE1bXHqrx0o+kY0tsvWWOybLplveAZk7bfRlsu/FOEPBzQrbZLuBXLRnquQG9myWDXo69fJ24sYHW9aG04lO4H5DdY5r15WsMLx1n/cm20056AVc2ONOuTbBW9ehHS8yE58fHecks8dc8BsBiNiGTQBUdzgSGecVa0VCVOTHv+t0H3eSKZN/5b1/EfPsi24PgTFkij59Hi6ydvObaOQ6o4nuOet7Ru1Zx/tPC2hMl+bpgAlz+Rxm9t56iXD3B8l/fIjzzXAr4Q2Ln55k+ttgu6VEIoQd+oPiFAN02mWIURFEVqoFBKoMdGUZnhGbOsXV6gWy8wNylmuLaFrK6Try1dqB98Snc99urVgQUIyBC/+KI4tYC1VJNXjaIwLTO8NFgRdmscm6YAf2sQUSpg6WfNRTWM5e1b+GjCgalzBuMi6hVQDB14PA/q39EWYtsjBBjiMfmqU+UNH1DdIJa8D1h81yfmIFEWHo2EK1gvGK8UtysCcMC+5YLxBcuPXRT2medZDmcXMHPdQ2lHuh2S0lAikBWeMqiwdmAqrTLcVcZ9ThDs4zQa/dZyNeFkOUM+gsU8yVZWcCVawcW8Cnc95kpi3ZGW9R24oM1sDgPITC4pFRvajfF9o0lBkFEEasYUUIwrLmARoN1ga3ck7vATdunsGFn31UjSlZ46qxtP0rVpHA/JJI5KHLC0gCMobg+pQxxZ8E4zdpx0TF3aG6QJqLWEDNDKAzqDCEz6Il5nJwnPPdCasEfIMkc2suJTqiWhPGpiJ0KpgEzNWgfji9scWZunYVsShUt05C1W2LWBeuLJaPNEtYyoP2lHjOL7/fo54ZcBF6+eiABn8J9n8XxuJ0a0rGLC+2Lf1CQjZT+MwWTUxlxLoB2E1lN25+3u29FvUWka617S921GgB8bLt3tIioMUhI0X6oVPHzBXbikSpg10eoadea0ejwiyVubYq5Om5HNhkDPqCFI8wV3dh4SzyzQJ49gT5/8aEYWvdIECH2u/kJOcR+RAvBjgwqIKIslhPO9NZZcBOCGhq1jIqCG/UAIws4ExllBY0t29ejtltiSswRv0imeiAt+BTuB8051ArqDLZR5i5HijVh84LDD2PXPwuY9oObGAGrmCy2i0tGAdN24agKMQoxGHzTdv/EHNo/lBwG7WaiVosZBZDf3GrrYSyEiDQBqSOrP7TI/Isl7uYEaXz7Na0wW1PCwgA7adDM4hdK3IVzyN+/mLpoDkIIqMhOV4wZGzRTQqFoHinKdgJhYTx9U+/smLaSbbKYta24a2ZIVGGjzIgTITrBD4TKG4zPkbCAixG9+sq+jpJK4X7QrN1Z09tWSr7mcZPA4FrO+psc0xXd6Y81DYS8XdY39gINEF3A2rZl3tSOGARtDDSm/Y9ogJDC/VCIYFaWiQsDNh6znP3mLaTx4APiQ9sl1xiyac3yqGJybsj45CKDSxPsqEadQaqAFrZdssLHNuCXB2RbK/jLLx/2Fc489R5TeyQqEsBOBW2EUCoIXDi2yjsWX9oJ9u1d0xq1LNgJy27Ei+Uy37AnmVYZzdQSpoqthWZOkGhBS0pZJhMhXL66b7+0U7gftBAgKHbqKW5E7LhGas9gVJNt9BidyZgcNzQDyDegngevII0j1oamjPgsokGgMaAgXrBjQ7YpZJsc+uSJR5XkOVrmaGZZuOiRcfe2u33L1dbeWrRXINOG/rfXmJ5fYPVtfdA+vVsRXwq2UiS2H76uv9nQDBXjL/Dmjy0Q//aZw73IGachYDen2GaAGwMi7RDIXkDKwFxWMWdvr/dkJTIwFaU0LOqYm2YOI5HCtHNQvhVP0FDi+4IbCSEX1FhiVpIPT1AWOfGFl/alBZ/C/YDpZIoZT9G6wcSITGu0qjF5Rrk+oriS45f6TI8XSIDJMcvWOSHmoFODWkMoFQFsLZga7ETo3VCKtUC56mFj667nkew90+9D4zFbU4qb3UtL5M5uMlVkawxFjg56FK+MWR7lrL61x9Zpy8pXJ9SLGVtnLOtPQLYJg8vt6At197uIa3LPVJGNLex0mWysxFyIueLmGvKioY6OG80Qmyl9U7Fs29daJp5SGkppGHRbZJamYS6reKZ3ktGopK4NzZJjctIABje29J44wanPOfRr39jzS0nhfsDiaARXFTM/RJuGsLHVLlFgpN3ZJctw13OGFwswht6pJdTMMT0mxKztzsk2BVHaYJ9CsR7JtyLFWkN2fURcWz/kq3w06WSCmfbRIie7uo6urrWrgpYl4tolYneCvqrbfvYihxhZflrZfLzH5X/cw1aQbyrLX4NsEtthko1iLl4ljZvZf1rVuInHjS3TZSGuNJxe3qB0nvlsyiRkhExYtGNKacglUIqnUUvfVGTi6Ut7O2ennC3XeGmyxEZdslkX1N4xmuZMNkpikTF6Yp7e12XP33GncD8EcTxGQ/sy3e5v09hNdKobZDIBukWpphXLcor1JwdUCwbRri+wBlu3+7K6ccRtNWQ3R3BjrZ0Rlxy4WFVI0yCjCfHqKzujIaSuMfPziDFo3bRBn3dr7qsS+zn1co6bKMeeaf9fREv7C7xp691/YYNw4+Z3+tHJHpGiQMqiHboqUC8oeb99PS3kE1aKLY5l7VpAm7HESsRKu1xIKe1ERStKZtt9GB7LVtkoSsZzBWuhz6V6mWvVPNemQ66U89xoFhgfc/Rdtud97yncD8n2HpuvEQO6PSY6BKgqzGjE8upJqgvLTFYyTOgCvopkmx67VWNXNwgvX0sjKg6TKrq5hY7vXBVQq4p461a736bGnU0dpN9Hihz30g3cK+16QGoN2i+IhUNCxEz9zi+LZP/Z4yvEY/NAu/lNs+JZ6bf94XV0rDc9fLQ80b9GUMNGKNkIJW/Kb5BLRd90Ad+t5jqQioZNNqOyFnPeVl7mql/gYr3C08UZ/nJ9QDOftePrU7jPkLu9Desej+MxPP8Cxa117Pecp1loN4Gw04DbrDCrm4QrVw9627bkdXynJXzvWH7C+3a7troBje0IjX4fMz9sd/oJEbMeiTdvEUbjNInpAGmZMz09167TVAp24CmcZ76Y4iSyXpdMbIYzK6xkW6y4TUrTcNUvsGkqjrPJMVORCeQi9E1Go4EgHmNqSg001rJq5gjdeMvoQPo92OPln1O4HxWqhJuruG8o5k1niKXDTpo22F+5noL9qFG9o6W2/UvBzA933gHEzc3DOrtHloSImnakUjMA68LOmk8AuQ0MbN19sDoHwJIbEbrF3S2KRRmaBlBMrCnFkkk7gaUhYlFK0zDvKkSUUIKeOg573O2Wwv2ICbduwfoGdtCHEPCTSRr6OCPieIw2HskcsfvcJTlgkynGt1tW+j40laMKlsIaoonUwQJtF1pQ4RUdMo0ZfVOzGUuuS+BFOcb3F5cYmprSBgrJKCQjEgmhHUmTSftuzGUB31PqE4M9D+MU7kdRDKlVN6O0qdPnJodIqxo7Cfiy/cBbJ5b1UY/CBgrnMaL4aBiR06PBR0tUQ+MslTrGoQ3+Umrekr/CUMYEnWBF6IslF+la9xEjEWsjwUC17Mic29N34CnckyRJOlpVSIjtXsQRxBum45xbrtfuneAUI5GoQhUchfVU0WG6WeGZBNZ8n83Yw6KM1ZJJYBqVqXjWomOqjoAhdl05CjQ9085eT+GeJEmy9zQEJLb94KFUpBFiY9jaKokqnBhu0XMN4DFWmYSM3HiMODIJLGQT+rbGEplqxrJMyUVYNJZMLIaKqVqmMSMiGBNRpzQDwfT7hD2cqZqmvCVJknREhFhYQgF+EIllu5ifMe3GHKujPlt1wdjnTH27/4KPlqxbcsBK7NadUWq1NGoIuz4TWzSORVMztBMGtmJQ1MRCiZkgg/6eXstdw11EzovIZ0XkaRH5uoj8Wnd8WUQ+JSLPdrdL3XERkd8VkedE5Ksi8o49PeNkT0x1zFf0L/i8/hmf1z/nRX12+yGb6np0pbo+INOupx9KRYuIGTSYLBBj28p2NlIFSxUc0+DwsVuN9XX2PhtrwVos2FSh0YjBkImlL8q8mbLkxszlNWrbjVp2Jrbt1aXcw3M88K9U9W3Au4BfFZG3AR8CPqOqTwKf6b4H+Bngye7rKeD39vSMkz0hCE/yQ/yE/HP+Af+MSzzPlm4AnCbV9chKdX1AMWJCu1cqBqyLZIUnyz0xmnbPnW57vZ39FNSw5XMatTtfAKNYsBb7rIaS1RhpNGAwDMRwxq1zLr/JwNXt8t5GCMtze3opdw13Vb2iqn/V3d8EngHOAu8HPtI97SPAz3b33w98VFtfABZF5PSennXywArpMd823nCS0WdIxQRgkVTXIyvV9cFJ1Hb5hyxibCTLAlkWKLKGzMadyUciSkSIKtTRMQoF05jtrPFuuzHtAaFWw1gbIpFCHEPxHHNbLBVjyCKhhHqxeM22nA/iu+pzF5E3AT8CfBE4qapXuoeuAie7+2eBl3b9sUvdseQhNdERm6yxwDKAS3WdDamu98EYojXEUnGFx7lImXnKrN3icqGYslRO6Gf1Tsu9jq7dnKMpabQdGhlVMBIxRErxWFEaVcaxYawNDUJG4LHeKq7nCTlg9nZb+3sOdxGZA/4Y+HXV9n3eNlVV+O627RSRp0TkyyLy5Yb9240keWNePV/l83wvb8fJnX1+qa5HV6rrfbKW0DPdB6m3/4m29yu2JpKbQN/VDLLb/w7TkLHVFIx8QaMWI7qz/G8pgVK03ZwDZapK0w2DnLNThnMTQk/brRj30D2Fu4hktMH+h6r6J93ha9tv37rb7ZWNLgPnd/3xc92xO6jqh1X1nar6znZDsuSgRY18lc9zisc4ITuNNZ/qerSlut4/cY7pokXKgDGK6ZYeMF24R70dwM5EStvQd3U3HDIyCvlOnzu0yxEYlFIE233oGoGAELv4PTYY43sQSsHs4Yeq9zJaRoDfB55R1d/e9dAngA909z8A/Omu47/cfQr/LmB919vB5CGhqjzNlxkw5IJ8z+6H1kh1PbJSXR+MLAzZfMwgLuLc7TEw2y13VdkZGWNEGbiaga0pbTt7tYmW8asCPiLUqjTdm6Wg0KhhFAuiGpaKMWEuEHKzp8Mh72US0z8Efgn4OxH5m+7YvwF+E/iYiHwQuAj8fPfYJ4H3Ac8BY+BX9uxskz2zzk2u8iJzLPAF/RQAT/ADAFeA96S6Hk2prg9m8sRxJufaVruI4mxARNsRMrRB7aPBiWC6sI4ImQlEbW8L47HEnZZ5QGgUQLFAg9Co7Y53vwSKSHSm3Ypxj9w13FX1/8DrDOJsvft1nq/Arz7geSX7bFFW+Gl+7rUPKEFVU12PqFTX+yfOsfF4jhZNO6bdRHIXsKL0sgaAwnpK22C6kTJNtN1yBG2QN9ESVQgYghpGmmO0HTVjuslMY3VsxJJxLBiHnDo4TBbwRYbWzZ5dT1p+IEmSBLCnTjI+KWAV59oJS4UNFNaT2UDebXpd2oasG+7YaLtGjFfD1GeMfLdiJIZGLaOYMzA1pal3hkVOY87NMMeNZsg45pSuYW44ZXy6h/TKPVvXPYV7kiQJEE4uEgoFb4hRiN1gmXaEjKe0nsJ6erahibYN9JAx9jkTnzGqczbGJS/eWOJvsrP08oZj/RF9167yaURZyKasFFtEFSYhZxIyogr9oubGgsLyItxc3ZPrSeGeJEkiwvREH4kCEYK3+GDx0dB0Swxsd8VMumGPG1XJ5rRgUmXUkwydWszI4iZCHAuTAJf9CqYB04B03TKhFCYnlWYxkC1U9Ho1deNQp8TFwZ5dUgr3JEkeeabXo563bQh7QwxCVTsUGNdZO0pGBR8MMRqCN/jGomOH1IKdGtxYsFOwFZgajFfEgwnd/QC2UUyjLH1TQSEWBesXBugyZAGiM9/xA87vVgr3JEkeeTKcw5eC8WAmQigsVWOpKNrpXlF2biUIeME0YCvBTgRbg6naFrptXhvq25vam0YxTcTWEVMF7EZFcaPk1vf1aAaCWknhniRJslek3yPmIAGyLYP4rA3Z7YmqCsYLEtj5Mp42xBuwtba/GLy2m3x0wb79RWzXrGm/b4PdTD1mbRP37Zc4cfMMq+84tntS7ANL4Z4kySNP+2XblVK33SNuLKiAREC5HdhNF9KBnda5BL2jlS5x9612x9sWu4SINAFpAmZrSrh+o9396dlvs2wMxMhejXRP4Z4kSeIDxUZEogGB6GiX/e1a7RJ3Bbnvulm61ji63ZrXNtC1ux8U8bE95iOm8hAUCQF8gKpGm3ZbPfWe8MxzSLZ3kZzCPUmS5MYq5Y0htnKEwqC2XWMdvd3d0oa83g7u7XCX7nYn7LuQb0Ib7j6CD0jjIUYIEbwnbo0g7mqnx4BWBzhDNUmSZNbFrRF2q24/0AyKWtkJ8+1+coJ2/e9teLd/8PYxtAv4JiDbId7tyUoIUDdobI9pVRO3tvb1mlK4J0nyyNO6xm5NcM5gaoMa6VrosQ3vGJGgt0N8F/FdkMPO4xJi20r3oQ3zGNtw9x6ahlhVr/l79loK9yRJElXY2MI4izHm9jFViNoG9fYxAJHbLfgYd74nxO62C3Tv2/t1g9Y1GsK+h/q2FO5JkiSArm9gsgycvTOAtwM8RjCmDXLZNRp9O/hjbAM9Kmgb7Nr4Nth9c2Chvi2Fe5IkCRCnU1i91S7eBbdb49D2mcdd4WwtWAOy3crf1WIPEZoGVUXr+sBDfVsK9yRJkk4cj2EyATFIt+2d7nx4GrpWu0GsRaxpQ35bF+YH2fXyRlK4J0mS7KYKGtD4Ro8FtOF294yY2633h0QK9yRJkvu1Hea6d+PT98o9bZCdJEmSHC0p3JMkSWZQCvckSZIZlMI9SZJkBqVwT5IkmUEp3JMkSWbQXYdCikgJfA4ouud/XFV/Q0QeB/4IOAZ8BfglVa1FpAA+CvwocBP4F6r6wj6df3Kfgga+wv8mElGUE5zlLfL9ALmIfJFU1yMp1TXZdi8t9wr4KVX9YeDtwHtF5F3AbwG/o6pPALeAD3bP/yBwqzv+O93zkoeMwfAO/gnvkvfw4/w0N7nKut4EOEeq65GV6ppsu2u4a2t74eGs+1Lgp4CPd8c/Avxsd//93fd0j79bRPZqz9dkj4gITto3btq18jpDUl2PrFTXZNs9zVAVEUv7Vu4J4D8CzwNrquq7p1wCznb3zwIvAaiqF5F12reCN/bwvJM9oKp8kU8zYYtzvIUecwAh1fVoS3VN4B7DXVUD8HYRWQT+K/DWB/3BIvIU8FT3bfVp/fjXHvTvPKJWOPwXkr3It95ykW+9DDz5IH9RqusdDru2+1XXrU/rx29y+P9vD8th13W3C9/pge9qbRlVXRORzwI/ASyKiOtaA+eAy93TLgPngUsi4oAF2g9qXv13fRj4MICIfFlV3/ndnMuseFiuXUT+HTAB/nWq6954GK5/P+ra/b2Hfm2H5ahc+1373EXkeNdiR0R6wHuAZ4DPAj/XPe0DwJ929z/RfU/3+P9SfYiWSkuAVNdZleqabLuXlvtp4CNdv7sBPqaq/11Engb+SET+PfDXwO93z/994D+LyHPAKvAL+3DeyYNLdZ1Nqa4JAPIw/JIWkae6t32PnFm+9lm+tnsxy9c/y9d2N0fl2h+KcE+SJEn2Vlp+IEmSZAYderiLyHtF5Jsi8pyIfOiwz2evich5EfmsiDwtIl8XkV/rji+LyKdE5Nnudqk7LiLyu92/x1dF5B2HewX3J9U11fUomqm6quqhfQGWdkLUm4Ec+FvgbYd5TvtwjaeBd3T3h8C3gLcB/wH4UHf8Q8BvdfffB/xPQIB3AV887GtIdU11TXU9enU97Jb7jwHPqerfq2pNuxDZ+w/5nPaUql5R1b/q7m/SDks7y53Tvl89Hfyj2voC7XyC0wd71g8s1bWV6nrEzFJdDzvcd6Y+d3ZPi545IvIm4EeALwInVfVK99BV4GR3fxb+TWbhGu5ZqutsOup1Pexwf2SIyBzwx8Cvq+rG7se0fX+Xhi0dQamus2kW6nrY4b499Xnb7mnRM0NEMtr/KH+oqn/SHb62/fatu32lOz4L/yazcA13lep6JK/hrmalrocd7l8CnhSRx0Ukp50d94lDPqc9JSJCOwvwGVX97V0P7Z72/erp4L/cfQr/LmB919vBoyLVtZXqesTMVF0P+xNd2k+bv0X7Kfy/Pezz2Yfr+0nat3BfBf6m+3of7bKqnwGeBT4NLHfPF24vq/x3wDsP+xpSXVNdU12PXl3TDNUkSZIZdNjdMkmSJMk+SOGeJEkyg1K4J0mSzKAU7kmSJDMohXuSJMkMSuGeJEkyg1K4J0mSzKAU7kmSJDPo/wPrtD3aVZqtLwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 3 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.subplot(131)\n",
"plt.imshow(np.sum(lul,axis=1))\n",
"plt.subplot(132)\n",
"plt.imshow(np.sum(fiss,axis=1))\n",
"plt.subplot(133)\n",
"plt.imshow(np.sum(lll,axis=1))"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [],
"source": [
"d5lof = binary_dilation(fiss,iterations=5,origin=0).astype(int)\n",
"d5lof[llung==0]=0\n",
"d6lof = binary_dilation(fiss,iterations=6,origin=0).astype(int)\n",
"d6lof[llung==0]=0\n",
"nonfiss = d6lof\n",
"nonfiss[d5lof==1]=0\n",
"nonfiss = (nonfiss > 0).astype(int)"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(17326, 13792)"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.where(fiss==1)[0].shape[0],np.where(nonfiss==1)[1].shape[0]"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"13792\n"
]
}
],
"source": [
"sample_size = np.min([np.where(fiss==1)[0].shape[0],np.where(nonfiss==1)[1].shape[0]])\n",
"print(sample_size)\n",
"ideal_size = 10000\n",
"sample_size = np.min([ideal_size,sample_size])"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7feac6a5d400>"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAACCCAYAAACpfRYBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS7ElEQVR4nO3dW4wkV33H8e//1K0v0z2X3Z3x2rter7ED2nAxxgGjEJIAVhxeyAOK4CGgCMkvRAIpD7ESKVGkPIQ8gIQURbJkFIhQCDJEgYQIX2LhXGxjB7DxJdhr7DV7md313Ptel38eqmY8u9je8e709EzX/yONpvvUzPSp+Xf/urrqVB1RVYwx5eRG3QFjzOhYABhTYhYAxpSYBYAxJWYBYEyJWQAYU2JDCQARuV1EfiYix0XkzmE8hhkNq+14ke0eByAiHvAccBtwEngM+KSqPrOtD2R2nNV2/AxjC+C9wHFV/bmqDoBvAB8bwuOYnWe1HTPDCIBrgF9sun+yaDN7n9V2zPijemARuQO4A8DDe0+N5qi6YjZZY+kVVT1wub9vdd2dXq+uwwiAU8DhTfcPFW0XUNW7gLsAmjKj75MPD6Er5s26X+858QaLL1lbq+vu9Hp1HcZHgMeAG0XkqIiEwCeA7wzhcczOs9qOmW3fAlDVRET+CPg+4AFfUdWnt/txzM6z2o6foewDUNXvAd8bxt82o2W1HS82EtCYErMAMKbELACMKTELAGNKzALAmBKzADCmxCwAjCkxCwBjSswCwJgSswAwpsQsAIwpMQsAY0rMAsCYErMAMKbELACMKTELAGNKzALAmBKzADCmxCwAjCkxCwBjSswCwJgSswAwpsQsAHYD5+EqFcQf2UxtpqQuGQAi8hUROSciT21qmxGR+0Tk+eL7dNEuIvLlYu74J0Xk5mF2fhxIEOLqNaQS4Wq1HQuBp/VxfqDf5WG9d6Mt1gHkM/9YXbeBq9XwpibBeaPuyuvayhbA3wO3X9R2J/CAqt4IPFDcB/hd4Mbi6w7g77anm2NKBFetIF5RBidIFOUhMOQnzdUc4d184IK2l/g/gDWr6zYQQRoTSKOBq9dG3ZvXdckAUNWHgMWLmj8GfLW4/VXg9za1f01zjwBTInJwm/o6dlwUgZML2iTwcRN1vIk6EoRDe+xpOUDAhX//PKcBFoq7VtcrIGGIBAFkGW6yCVLUWeSNf3GHXe4+gDlVPVPcngfmits2f/yb4S7694sDz4MgBCe4ehWJoh3rzoA+QFzctbpeAe33QQStRuB7+EeP4F81h3/V3I7W9FKueCegqiqgb/b3ROQOEXlcRB6P8yde+WTZBXelXkOaDSQK8zCg2EoYwbuG1fUKiUBQ7M/pDyBN83D3fbzpqZF2bbPLDYCz65uAxfdzRfsl549fp6p3qeotqnpLwO5JxB0jglSijRc6vg9hAKrge/kygDDA1XbmM2SY1yHIu2d1vRIuikAVaXfB90lnmsTXzZLOTaHTzV1zxOdyA+A7wKeL258G/mVT+6eKvca3AiubPiqYTVythlQq4PtIvY40G/mLH/J3j/UnSJoh1cqO9OkAVwPsK+5aXa+ATDbRMECnm6Szk6QTIf7ZFbyFNZLp2q75GHDJGBKRfwR+C9gvIieBvwD+GvimiHwGOAH8fvHj3wM+ChwHOsAfDqHPe5MIXqOBNBtovQpOSGshg+kKXj/D6yW4dh+SNN90jBPodiFJIMkQ30eTZNu681N9lCXOE9PnP/XfuJ5jHOGtnOC5pog8j9X1skkUwb4pslqIJBlpNUDSDElStFbBX+rAwVl44aVXQ39ELhkAqvrJ11n04df4WQU+e6WdGkcuivIXf62CBj69Q/k7vt9Lcf0E9YSsGuBWE6TbR30PpidhtQWD+NIP8Ca9Q9732guU51T1lguarK5b5mo1uP5assAjmQhRJ6SRo/byKjhHFgW4dg88D1etknU6I+3v7vggMu6ch5s7gE7U0MAjmYwIF3t47QH0B2gtIqsEZJWAwXQFv52QVn0GUz7V+QbBqUW8WpXkzFnI0lGvjdlEfB83PQ1TDbLJGu1ranRmPbr7Bb8HzZdSame6pI2IZG4CN0hxgyQ/PDgzbQFQBt70JNnUBBr6ZL4jeKWDtLtkEzUk8JFUwQlZ4Aha+bt9dGaV6JTSPzRJcuMsfjdlcNM1VE538OYXyFZWR/7kMeCmJtG5faTNiLTq057z6M8Ig2ml/gxIpiT1gOj0Kt5Kl+6RKfCE8IVzaLOOq9fJ2u2R9d8CYAdou4O8PJ8P8Al8NAqID07nm/2RR3/Kx+9mBGsJSS0grTpCz+G/0iI6vUraqBBPRYRLA9KJEJmbwUUhsrxCurwy6tUrt9l9DObqJDUPyZTO1UK4DFc/lFI920U9YTAZ0n7LNC7OWPjVgPYhj9nHDhO0M8LZCfyHn87HDYyABcAOyHo9nOchSQSeA+fw2n3SRoXObIDfV7JQiJs+LlaSioPJgMHUNHHNEa5lVM60IFVENd9x5BxSreIGsW0JjJB0+7hBHQkdWSgc/J8B0enWq8uBKFXSio/6wsH/btG5usr5d3lUzzta7/fwP/geDj4yoPLiIiytkC4u79hHPQuAHeA1mzC3H3UuH/obJyT7m7QPRkyc7uO1Bmjg5U8SJ1QWY1avjZAM6mdj0tCBKq7VyXcOeh7qBOpV2xIYseTFE7gXT2yMeHCVCnL0MBoFaOCR+S5PAVWSqs+g6VM92+fw/cLyDRXCFcGlcOo3Aqpvu4rK4iwzP14iffb4joSABcAOSNfW8JsNspkGWSVAPUdwrsX0C+fyQ4K+B60+nufQ0Ed9x8xTfXpX1WjPBVSWUxbfNY1Lp5h8vo23sJYHAKCBD/tn8DyPdOHiUzbMjvM8ZLWNOEe2r8lgtkJcc2S+oB4kVWHtUI0sAMlg+vmUzBNcorhEiRYGpI0I7y1HyF46icaDoXbXAmAnqJKcPJWPoCd/Q0hF8Pbvh3o1P+bv52f/SZyiIogqtRdXqJwN6R2oIqq4BBbeMUHYqhOspUQLPfzzq/nvV6u4RoNsbW1062mQwIcoJJ2ug3NE53tU0gyKUd8aOCTNyEIfDR0SZ6jnyAJHGjnUd6SBQ2eb+HKY9PhLQ90SsAAYtSxD4oS0XmGwrwJOUAH1pXjXEDIPwlaG38movKJ4/RQ3SEGVrFHFrbQhzXATdbTXH/q7hrkEVZJmhNdNkH6Kt9LOt9gCH818kqkK/nIPN99BK2F+UliSopFPOhEVYwc8sqsnCYMb0BdOkPV6Q+mqBcAISBDijh7eePGr75FOhIRLPVxn8OrJP4MYwvyzJKnSu2aCwbSPSzwQkJT8EOKhBvGER1IRwtYRGg8dJ31l4Y07YYZCixF//an8TIhwoZWf9OW8fFh3nCKDjMV3TtF8uYK/0EXiJP/q9XGtHulkHa8bF2NGKvhHDiE/f3kowW6XBBsBN9mgf3hqU4OjdyDEW2oj3T7S6eWnklYipNND2j0kTanMt6md7uH1lep8D6+f0Zv2WLkuQB00ftGnct7e/UdGBLd/hmxfk9VrPcKzLSROIEmR3gDpD5Bun+DsCjNPLBHXfVaPTZHsq5M1qqT7m2S1Chp5+cliSQYKyUwdb3b/ULpsWwA7TQRUCVb6SKYgQtasUjnXv+D0YOn00Mk68TUz+IttiBNEBG+tDwKrR2uoAy+GiTMpkimZJ4Sr/fz8AbPjJAzRSogGHpMnkjzIIT9sm2UbpwRrNUJ6MbUXl+kdnmTxWA20RnUpI6kIXl+RLN9huHK9I24oLjnC9d+cJHvi2W3tswXATlPNnwjFSSAaBriVDm4FSLP8MGEREm65hXRDkn113CBFnZA0QrLAMXG6T+Y70orLdxpmShY6JMnIusP5vGjemKvVIE5wrR7RQvHSErnwug+qSKsDUYjWq0TnOsy0QxbfVqV10GP/k10GUwGtqz1WboBgDeqnBBTU3/4NdguAEdA4wZtfQgdxfh1A59As27iPCIhDfA9RJZhPiK+aJKn5+N0U10+JmwFQ7AfQjO4Bn37TsXZoiql9b8d/8EcjP9OsbLTbxfVqaBQSzK+gi8uoKlKpIH5+MZCNMOgP8s/9UQhZxswzytrRKqc+WMXrQ7imzDwFQTfLDxHGijsxz3YfD7AAGAXn0PXx3+JgqglrLdCMzRXWGIhjpFLBP7eKN1GlP1vDxRnVk2ukExG9AxGd/T5eDI3TCW6geN0k/7tqJw7tpKzfR+I4P89j/tzGnnsZDHDNZh70gzgPgzAPcFTJaiGDmRC/q+x7Nq9Z5oEouFiRFGovrQ5lx64FwCjEMXivXvFHWx3ShcX8IiHBRSVJErTVgq6Ha3epxCn9qyZYOTaFixVvoDROJfjdBK+dn0jknTxPYmcN7jxVdK2FdroXHLbTfp9saakI5WzjYiBSqyFRiP+LV/DP5RdoVc/lZ4dGPpJmuF6yESjDYAEwAlm/n58bEPjoyipZv5+/E7Tb+anDYZBfP865Vy8ZnimEAVno4/op9TMZ6gR1guunGy9+1xvYIcARer3zMjZfzEWTJN/LP4hBMzRJcLUartlAPC9/4a9kZAtLpO2ODQQaO5te7L9U3Cwl673aJr6PhGG+X2BlFbeyilepwPo1A1XzqwiJoHFMtri8rVcOMkOiesFx/fXgcMWl4XSttSOjOi0ARmkLya5JcuELWgQXJ0i/ePIUWwg7eQaZGY6s00HjBAl8sm53Rx7TAmCvUSXrdHCQXzlYlXRxyfb4jwmNBzs6lNsCYC9aD4E03dh/YMzlsADYq1SHdoKIKQ87F8CYEtvK9OCHReRBEXlGRJ4Wkc8V7TZF+B7W0w7/qz/gYf0+D+u9vKzPry/yrK7lsZUtgAT4Y1U9BtwKfFZEjmFThO9pgnAj7+T98jv8Gr/NSV6gpasAB7G6lsZWpgc/o6o/Km6vAc+SzwxrU4TvYZFUaeZv7vgSUKNBny7AFFbX0nhT+wBE5Drg3cCj2BThY6OrbdZYZpIZAN/qWh5bDgARmQC+BXxeNd9WXHc5U0nbNNK7Q6IJT/Iwb+UmfAkuWGZ1HX9bCgARCchf/F9X1W8XzVc0RbhNIz16mWY8ycNcxbXMysabeWJ1LY+tHAUQ4G7gWVX94qZFNkX4HqaqPMPj1GlwRH5l86JlrK6lsZWBQL8O/AHwUxH5SdH2p9gU4XvaCgvM8zITTPKI3gfADbwd4Axwm9W1HLYyPfh/kV/K/rXYFOF71JTs5yN8/JcXKKmqWl1LwkYCGlNiFgDGlJgFgDElZgFgTIlZABhTYhYAxpSYBYAxJWYBYEyJWQAYU2IWAMaUmAWAMSVmAWBMiVkAGFNiFgDGlJgFgDElZgFgTIlZABhTYhYAxpSYBYAxJWYBYEyJWQAYU2IWAMaUmAWAMSW2lZmBKiLyQxF5QkSeFpG/LNqPisijxXzx/yQiYdEeFfePF8uvG/I6mMuQasoP9QEe0ft4WO/lBX16fVFodS2PrWwB9IEPqeq7gJuA24upob4AfElVbwCWgM8UP/8ZYKlo/1Lxc2aXcThu5je5VW7jfXyEBeZZ0QXI5/yzupbEJQOgmA++VdwNii8FPgTcU7RfPI/8+vzy9wAfLuYXNLuIiOBLPjGUkqGvTgLcwOpaGludHdgr5gU8B9wHvAAsq2pS/MjmueI35pEvlq8A+7axz2abqCqP6H08xHeZYZYqEwCp1bU8tjI5KKqaAjeJyBTwz8DbrvSBReQO4I7ibv9+veepK/2be9R+4JUR98E7wXNvOcFzp4Ebr+QPWV0vsBtqu+7IazVuKQDWqeqyiDwIvB+YEhG/eDfYPFf8+jzyJ0XEByaBhdf4W3cBdwGIyOOqesub6cu42C3rLiJ/DnSBP7G6bo+9sP5bOQpwoHjnR0SqwG3As8CDsDG97MXzyK/PL/9x4D+KmWXNLmJ1NbC1LYCDwFdFxCMPjG+q6r+KyDPAN0Tkr4AfA3cXP3838A8ichxYBD4xhH6bK2d1NchuCHERuaPYdCydcV73cV63rdgL678rAsAYMxo2FNiYEht5AIjI7SLys2KI6Z2j7s92E5HDIvKgiDxTDKX+XNE+IyL3icjzxffpol1E5MvF/+NJEbl5tGtweayue6SuqjqyL8AjH1R0PRACTwDHRtmnIazjQeDm4nYDeA44BvwNcGfRfifwheL2R4F/BwS4FXh01OtgdR3fuo56C+C9wHFV/bmqDoBvkA85HRuqekZVf1TcXiM/1HYNFw6tvXjI7dc09wj5eIuDO9vrK2Z1ze36uo46ADaGlxY2Dz0dO8UZdO8GHgXmVPVMsWgemCtuj8P/ZBzWYcv2cl1HHQClISITwLeAz6vq6uZlmm8j2uGYPWiv13XUAbA+vHTd5qGnY0NEAvInyddV9dtF89n1TcDi+7mifRz+J+OwDpc0DnUddQA8BtxYXFwkJB9d9p0R92lbFafM3g08q6pf3LRo89Dai4fcfqrYa3wrsLJpk3KvsLrmdn9dR70Xknzv6HPke43/bNT9GcL6fYB8M/BJ4CfF10fJT6V9AHgeuB+YKX5egL8t/h8/BW4Z9TpYXce3rjYS0JgSG/VHAGPMCFkAGFNiFgDGlJgFgDElZgFgTIlZABhTYhYAxpSYBYAxJfb/a5KpiXeyXKUAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.subplot(131)\n",
"plt.imshow(np.sum(nonfiss,axis=1))\n",
"plt.subplot(132)\n",
"plt.imshow(np.sum(fiss,axis=1))\n"
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.colorbar.Colorbar at 0x7feac69a74f0>"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAD8CAYAAADqmhgGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArqUlEQVR4nO3deZgcd3ng8e9bfc19j2Tdt2RLtuRDyMbYHDbBxmQxLITYSziyziPCkUDC7rOGPA+wSXgekieYXZbDUdYOJmswDpiggME2xgSM8SELWadljWzJkqxj7kNzdXe9+0fVSD2jmememe6p6un38zz1THVVddU7PV3v/H71+9WvRFUxxphS4AQdgDHGzBZLeMaYkmEJzxhTMizhGWNKhiU8Y0zJsIRnjCkZlvCMMaEjImUi8qyIvCAi+0Tkf46zTUJEviciLSLyjIgsz7bfgiU8EblZRA76wdxZqOMYY+akIeAGVd0EXA7cLCLXjNnmDqBTVVcDXwH+LttOC5LwRCQCfB14O7AeuF1E1hfiWMaYuUc9ff7LmD+NvUviVuA+f/77wI0iIpPtN5rXKM/bArSo6ssAIvKAH9z+8TaOS0LLqCxQKMaYXPTS2aaqzdN9/01vqdT2jnRO2z6/e2gfMJixaJuqbsvcxi84PQ+sBr6uqs+M2c0i4BiAqqZEpBtoBNomOm6hEt65QHzHgaszNxCRrcBWgDIquFpuLFAoxphc/Fy/f3Qm72/rSPPMI4tz2ja24PCgqm6ebBtVTQOXi0gd8EMRuVRV984kxsAaLVR1m6puVtXNMRJBhWGMyRslrW5O05T2qtoFPAHcPGbVCWAJgIhEgVqgfbJ9FSrhnQvEt9hfZoyZoxRw0ZymbESk2S/ZISLlwO8BL47ZbDvwIX/+vcAvNMtoKIWq0j4HrBGRFXiJ7jbgvxToWMaYkHCZWultEguA+/zreA7woKr+WET+GtihqtuBe4B/EZEWoAMvz0yqIAnPv4D4CeARIALcq6r7CnEsY0w4KEpyitXVCfeluhu4Ypzln8uYHwT+YCr7LVQJD1V9GHi4UPs3xoSLAukcqqtBKljCM8aUnlyuzwXJEp4xJi8USId8BHVLeMaYvMlbk0WBWMIzxuSFonYNzxhTGlQhGe58ZwnPGJMvQppJ790PnCU8Y0xeKOBaCc8YUyqshGeMKQlex2NLeMaYEqBAUsP91AhLeMaYvFCEdMgfk2MJzxiTN65aldYYUwLsGp4xpoQIabuGZ4wpBd6Ix5bwjJmURKOkX3/Z6GWuIk+94N2vZIqCqjCskaDDmJQlPBMMEc6+ZwtuVFCB4erR135EIbb0akSh7mcHSHd1BxSomQrXruEZc57E4jhrV6AiDNY6TFQgyEyC6bVLib7WQeq4PQcqzLxGC6vSGnPO8Jsuo3tlfErvabu8itiaSsrbFpD45R50aKhA0ZmZsUYLU+Kcykr6brqUoRrvRHBj09tPslJIVsaI3H4lAE3PtsNrp62qGyLF0GgR7uhM0XM3rKS/ySEdh3Tcq6rOxMh+Tl/XSP8b1uYnSJM3aZWcpqBYCc8UhESjRJoaaV9bWbBj9DdHqbhkDekDhwp2DJM7RUhquFNKuKMzxceJEFm5lN7L5jHQ5My4RDeZdBzOXNdEcySCnG4n3dpauIOZrKzRwpQc99rLaFtfPmvHU4Ez19aT6K6j7hdY0guQEmx1NRfhTsemqDhlZfQtLQvk2EO1gs5vDOTY5jwXJ6cpKFbCMzMngnPZOs6urGG4Krj/8K3X1NN0MI4mhwOLoZSpEvpuKeGOzoReZO0qBt/xOlqvrufsvIBvKxJIX7MBJNzVqrnKa7SI5DRlIyJLROQJEdkvIvtE5JPjbPNmEekWkV3+9Lls+51RCU9EjgC9QBpIqepmEWkAvgcsB44A71PVzpkcx4RTpKmRM9fPm/BuidmmAp0XlyGrr6H54cOkT58JOqSSk8dGixTwaVXdKSLVwPMi8piq7h+z3a9V9fdz3Wk+onuLql6uqpv913cCj6vqGuBx/7WZYyJrVzG0cXlokt0Ilel3bjYzowiu5jZl3ZfqSVXd6c/3AgeARTONsRBV2luB+/z5+4B3FeAYJiCSSBBZu4q2a+fRtXpqt4jNpp7rV+BUVIATsow8x6VxcpqmQkSWA1cAz4yz+vUi8oKI/FRENmTb10wbLRR4VEQU+EdV3QbMV9WT/vpTwPzx3igiW4GtAGVUzDAMU2hORQWyZCEdm5tIVUhB+9flw0CDw+D7N1HXMkT8tR7o6rEqboF5z6XNOZk1iciOjNfb/PwxiohUAT8APqWqPWNW7wSWqWqfiNwC/BuwZrKDzjThXaeqJ0RkHvCYiLyYuVJV1U+GF/B/uW0ANdJgg56FnNPYwKk3NQcdxpSoQOeaBKxppuFANY4lvAKTqQzx3pZxGWz8vYnE8JLd/ar60Nj1mQlQVR8WkW+ISJOqtk20zxlVaVX1hP/zDPBDYAtwWkQW+AEvAOxbVsQkkSDS3Ez7mxYHHcqMdK0pI7J+LZGamqBDmbO8xzTmrZVWgHuAA6p61wTbXORvh4hswctn7ZPtd9olPBGpBBxV7fXn3wb8NbAd+BDwJf/nj6Z7DBMsSSQYvHEjvYuLv7umG/UGHKg8XUf5j54NOpw5SVWmUqXN5g3AB4A9IrLLX/ZZYKl3LL0beC/wURFJAQPAbaqTD5E9k2/yfOCHfoKNAt9R1Z+JyHPAgyJyB3AUeN8MjmECNFeSXabB+gjcuoWqJ14k3TP2kpCZqXx1PFbVJ2Hy+rGqfg342lT2O+1vs6q+DGwaZ3k7cON092vCIdLYwEDT3Ep24A04cHZ+hKq49V3JN288vHC3ZtmdFmZcQ5tWkArmtthZkdywNOgQ5iBvxONcpqDMvX/hZsaiy5bQuTS8fezyoWdZGfWJhA0Xn0detxQr4Zkio/EY6bmd70iVwfD1lwYdxpySz3tpC8USnhlFYnF6Ly2u/nbT1T8/RnTRwqDDmFNseChTVPpvuZz+5tL4P5isFLTK7vLJF294KKvSmiIRXbyI4erS+koMrKhHovZ/P1/yNXhAodhf2pyTXNpEsiLc/6HzrXdJlLJoFE2lgg6l6HmjpYT7H6YlPAOAc/l6WjeUZvVOViwhcqaddHtH0KEUNe/WsnAnvHBHZ2aFs/FiWq+qDf0IKIWgAqevb6L7hjVE6mqDDqfIeSW8XKagWMIrcZF1q2nbXB+6gTxn22C9Q/t/Wk/qhquCDqWouUhOU1As4ZUwicUZWF6Haxc2AEiVw0BzDGfTJUhsjndELICRVtpcpqDYV71EORsv5uyKGs7OL/Gi3RhDtcKZq+uouWgj8Uef985ik7OwN1qEOzqTd5GmRuSKDbRuqbdkN4nepTGG33YVTtkcvqE4z/L5TItCsRJeiem7blXwj1MsAirQsyxG+eIF0PJK0OEUBQVSIS/hWcIrARKLE2msp+v65SXXsXjGbBipKQl7ldYS3lzmRIisWUHv+kYGmpyS7HYyU2eubaTpUBxNDgcdSvgFXF3NhSW8OSqydhUDK+rpWWYlFDM7imEAUEt4c5BEowwtqbNkZ2Zd2Et44a5wm2lJv+EyulZZP7K8EHDWrgg6iqIwMgCotdKa2SGCXruJzrXWlSJfVKB1SwPN7hrSBw4FHU6oKULKDXcZyhLeHOJsvJjW9eXWOJFnbgzS1fZPJBdhv4YX7nRscuaUlXF2RbUluwLpX1yBU10ddBjhplalNbOk75ZN9DfZ/69C6W9yGHrPpQw2CAu/+qyNnzcOe4iPmRXRRQsZqrE/ZaElq4TeS5Ik37SJyJqVQYcTSlbCMwU3uG7BnH/KWJh0rYoTnzePmkMvBx1KqChCOuSNFuGOzmQliQTpMvszzrZklZC64SobRmoMGw/PFJSzdBG9i62gPhv6FyiJk95nrQKvvi1O+ur1AUcVHloEjRZZE56I3CsiZ0Rkb8ayBhF5TEQO+T/r/eUiIl8VkRYR2S0iVxYyeGNmU7RPSFUqIwWURJcQ7R0KNqiQUZWcpmxEZImIPCEi+0Vkn4h8cpxtppxvcinhfQu4ecyyO4HHVXUN8Lj/GuDtwBp/2gp8M4f9m+kSIdVsXSVmS1m7Ej0r9C31Xic6FDcRBQl3y+Tsyet4eCng06q6HrgG+LiIjC1OTznfZE14qvorYOzjnG4F7vPn7wPelbH82+p5GqgTkQXZjmGmR+JxOkr0SWNBiXd5pTw35p207ZdVIREbX3BEvkp4qnpSVXf6873AAWDRmM2mnG+mew1vvqqe9OdPAfP9+UXAsYztjo8TJAAislVEdojIjiRWLZgOp6oy6BBKTqxPqT7s0LUp6S1woOsPN+NcenGwgYWAKqRdyWkCmkbOf3/aOtF+RWQ5cAXwzJhVOeebETO+2q2qKiJTHvhfVbcB2wBqpMEeHDANPW9ZY3dWBCAypEhZGoiiAslKQcus4QimdGtZm6puzraRiFQBPwA+pao9M4kNpl/COz1SdPR/nvGXnwCWZGy32F9m8iyybjXDVdbIHhY9q6uI1NQEHUaglPxVaQFEJIaX7O5X1YfG2WTK+Wa6Z8x24EP+/IeAH2Us/6DfenIN0J1R9TV5lK6vmNOdjRM9Snm7S3m7i6SDjia77pUOUlvaCS+fjRYiIsA9wAFVvWuCzaacb7KWw0Xku8Cb8ercx4HPA18CHhSRO4CjwPv8zR8GbgFagH7gj7Pt30ydU11N34LyoMMoqPrfniB11Ls8U73xYgYWV4emv2H/QiH6WuLc64Fm4dN/9BDf++VNyLHjAUYWvDw+1fINwAeAPSKyy1/2WWCpdxy9m2nkm6zfIFW9fYJVN46zrQIfz7ZPMzNOVSX9zcVbnZU0ONlKbe75M8fd/SIVyTX0LWwi6GfE9F8k6Ppeyn9TjVeJg6EmlztqT/G9YEMLhVyrq9n3o0/C5BcEp5NvwvEv08xp4kL18fOji5Qd78Xd++Kk7xk7Fkn6wCGql9XRszS4r2z/RV6XlJrHq84VZdSBxrXtgcUUJl4rbbj/EVvCKzLOpks4u7gq6DByEjurNOzqBNcdNVqwG2BM09W/UEgnlJrDo5cnb+rmmcsf4FMnryZ6pociuNxYUHms0haEJbwiM7Cwit5FxfFniwwr6X0Hgw5jlGS1V0tyY9B/ySCxYwkSnYIz7HU3GaEO9C2D1CKvj2j0tQRVR87vZ2Ce8K9/8mWaIy6b7v7vxLthYVe4ftcg5KtKWyjFceYYACLNzQzXFEevflGoPtJPGP7h964A9T+2dLlXBpOUUP+bBH1LoW95GkkJVUccogNexH3LwBkW6p9MjLtPjcCGeDl97iDxrtHJslQpuXc5CYolvCKiC5sYqg33F2pEw95+dMfe7BsWkkDPKoj2C5ICFGoOO2Rm4apXARxSlULPxSlIeBXusiNxKk5dmMTUgdTNXfzFuv8gqWmu3fHHRJPedmevWUXZj0v7el7Y074lvCIRWb2CM6+rCzqMnEXb+0jn8YKOU1lJuiz3ZK8O9C2HWK9Qfjp7HNGzSv0LERCvKCjp8d/TeWWKl173/4hJhHceegfRn9Uh/kXJgaYIlZesASD9Ykv4L2jlm4K64f6HbAmvSHS+bn7gXTKCJMsXc3Ze7tX53hUQPTtxskuVC/PeeYxXn15M1VH/GBO0pgw2CYtvfBWAtzceJeYnxZTrjHpPskI4uLWRyuMOi46ewO3vzzneucKqtMbMkFNRgRvP7auqDvQu90t2Z7xkN1wrqONdd/v8h++nxhnEwWVJtJvWFRX0uwn+6sCtuD9twkmd727S+8YB/s+W71Ln9LM6Nkhr2juZDwzDzqElHPvpcuIo6oAbF3rWpal+KUJ0QBl4ywYSP3muIJ9HmIW9UGsJz4STCJGLV4MI7Vc1ksrxxpK+5RAdOJ/sABa9/SgfWPQ0769u58BwP39/6iZaB6s4+d3lXP0nv+Mbi57m5qse5I3xd3P8dD0AldWDHLr6O/4eImz53YdJ/6hp1LHieIOB9i6HVF2KmgPRc40eycoI4zd3zF0j99KGmSU8E07icObapnOtq9kMzBMGliQpPx4bXY0VqIgOUyZJ+t1h3vHjv6B+9/lrA4/+8gpWNF7KOy7bw68u+yFcduG+Hx+I0PN8E+MNxtW9GuLdQs3LEcJ/yb7AFMI+fI8lPFPUUuVCsgqGa5XavbFz3UP6FwhuXHnjW3dz9+JfcyTVz6U//DT1e0ZfCK09BByK8ZudV3HHH8W4pWH3qPVJjfC3/3w7lR2jk9lI1TneJZS3lniiy2BVWpMXNa/007bRBvwcoQ70rFY0ojhJ8RKXX8LquljZ/p/vYkO8nM+3buAte9/D8ZZ5FyS7TM6wsuvey9g1ThEvwfm+eemy82d0rNuxZDeKWCutyY/I/iOwcUPQYYSCOtB1iVJ+2iHWB5HB80mnZyX84N3/m5i4rHrgT6k46RDvVur99YONwsA67+4JpzVObQvZa6ICvcvASUGizUua4kKi05LdBUL+kVjCM0XFjQm9K10qjznEe86fXYNNQv+yJM/c8r/YM1zD7d/6CHUnlcwzcLBRGGx2qXs6fu49PSuVaL9XKilr90p6Y/Ut8ZJd5XEI/RkdJLVGC1Oizq5rJLawlsSrHaRePpK3/fasdik/PTrZJauFP33/T/iz+qM8PVjGR/91K9UnRyemgWYhWavUHZRzOctryRWSVd6C7rUujNNIEm+LUHnCEl1OQv4xWcIrEuneXpp29dF2eXGMlNK3IAILIjRoPc6RY+DObBwRNyr0rPVKdrFerzvIYKOw5K1H+dKKh1gcTbHy0T+n4sUE1R3e+lSF0LfURRuHcc7EqTksF5yQ5WeUcv8BBanyyLhjgEfPhvwsDhUr4Zl8UMUZTAYdxZR1rCujMXkp8tQL096HRqB3lUv5ST/ZAR1XJ3nppn/knu6l/O3xd/D8oeU0PB1jJKOdXQRDTSkqjkUpO5Bbj7iRPnRmBkI+9pclvGLieqOQhPwyyQU611XQ8HRk/FLeBA+xPvesV4HOjWkQSHSe3/Yzr3+YB/vm8dX7b6W8VWnIeG/fEu99DbuC6xtXbH+jvLB+eCaf3H0Hqa/bROcl5WH/Xo3ixmDwHVdRtXP08x7c+hraXlc//pv8OxiSS4ap2psgNqZa+cuOdbyhvoXImELv2UUCjvqjoARnuEbgmo3w9O7sG88h1g/P5I8qzpO7qI9cQce6sqCjmZK+BRH63rEs9+2XAS7UPxVnvFLa88eX8MH5TzHUoMS7/PcswUt2R/MR8cyogEackF/RKgBLeCbfIr/Zw0UHG+l804qcxscThfmPnYDBoZz2n17cTOuV1TMNc8rcqJCqgv6NA8ReLvO7gZxX+Z5TfGT5rwC4qeJV7mp7PdUvn18/vCBJ3c7xE2QQuleV09gyj/TpM9k3nitCXvWwhFeENJUideo0lScWMFwzcfW26mQaxx+cMn3yNDqUW8KLilDT5N2tP1ztMFhf+HGpUhVC7woXrUgTPzI62XVudNGyNB9b+lveX+0NsDmkcX766vppP1h5NqTKQOJz+OHB45Bw/K+ZkCW8IuY8tYd5natREdo315OOQ8PBQaLtAwDooVfOJbmpfA9TJ08RP3kKgPKmRmouaiZdW0b7hsI8C9eNesmu6qhDrO/C7H33zffytorRF+u63WHkZ+ev/w02CdGOGOKG/Iyby1TAbi0zBeOmzz0kp/GVSnAc3LP9uDPs85Yp3dYObe2IE2He7gokEafrrWtIlgtuLA/7Twg9F6epPRC54C6HdEK46N1HeVN5PzD5wZI1SvSshC7htb9xMbXfORH+q/n5EvJf0xLeHOGePVvgA6Rxe3uhF6ofaCf51qvoXjmz6lqqXOhb7lL90oXJDqBnTZpdF/+EbMkuzNIJAXFAS+QBjpbwzFxU9sKrlLeUM7Siia6ViZzHrRuhEehdk6b6cGRUh183Kuc662t0/LOnzx3kjb/9UyozO7l6d4mZoFnCM3NRurUVWiFy9BhNT8UZeNsmNAIDDRHcLN+qVLnQu9qvxibPD8Pev8jlh+/6Csv8RBcjAowuRf6gr4a//erHqEzqqJOr+gh0rXNxhhzK2kN+1s1VRdDxOGsjl4jcKyJnRGRvxrIviMgJEdnlT7dkrPuMiLSIyEERualQgZuQUEWHhij792cp/7dnaXq2k+pjqUnfMrBASbSdT3YAfZcMc/h9d7MxXkatU06tU06Fc2GV+TM73+VVf8fmNPWGdk9VhCvZqQPO+jVBhzFrRHObsu5nnLwzZv2bRaQ7Iwd9Lpf4cinhfQv4GvDtMcu/oqr/MCaI9cBtwAZgIfBzEVmrWioXMIy790UqkqsZaG5GBdLTfLDD8VQfDU6co6nzydNFKH+miguynUDfYpA0VJ6YfuyFoA70L6+hLOBH9M6a/P2/+Rbj551Mv1bV35/KTrMmPFX9lYgsz3F/twIPqOoQ8IqItABbgN9OJShT3NIHW6g72EKkrpbBzavRiNC7NDal2s71D/8l77/mt/z0G9d51+b8Eyk6zhnVtxgQqDyWl/DNDOSrH94U807OZtJv8xMistsveo50iFoEZH7tjvvLLiAiW0Vkh4jsSJJbh1hTXNJd3cR+/jzxR59n3q9aqTydPn9CZJ4YwoVnioI7kiEnOYn6lnlvDfreWeNTyW2CppHz35+2TuNorxeRF0TkpyKS03Dg0220+CbwN3hfxb8Bvgz816nsQFW3AdsAaqQhXBdeTH6pkj7YQsXLcapWLKFjUxNVR87/r33/nz3CH9XugYzngklS+Ml9111QohuuPd+KO1yt4IavGluyRg8wnU2bqm6ewdF2AstUtc9vQ/g3IOvF0mklPFU9PTIvIv8E/Nh/eQJYkrHpYn+ZMWhyGFo7uOjJZqpe7Qdg4KIyFsfbmRe58AFF0f7RZ0//AmG43j1XNY71SuhLduJC5SvdlMxF7FkquqhqT8b8wyLyDRFpUtW2yd43rYQnIgtU9aT/8t3AyCXZ7cB3ROQuvEaLNcCz0zmGmZvSnZ1UP/D0uddtn7+WG8pfo83PCDuH6vjYj+6g+rjDyNmjEW+I9nSZUnfgwlGLw0xcSB9oCTqMWSOzNACoiFwEnFZVFZEteJfn2rO9L2vCE5HvAm/Gq3MfBz4PvFlELsf76h0BPgKgqvtE5EFgP5ACPm4ttGYyiU64dd8H6N9+EQCiSm0aRrLacI3QtzpFtDtCzcsUVbIrSXn6+0yQd2IAqno38F7goyKSAgaA21Sz37+XSyvt7eMsvmeS7b8IfDHbfo0BmP/Vp3C2lTH0gfnjtuKeXZamqiU66qE9Jpxy7WOXiwnyTub6r+F1W5kSu9PCBE5ViQx5wymZIlfsd1oYU2g6NETToy8T77VSXNHTHKeAWMIzoZA6dZqGJ14hMowNAlDE8nVrWaFYwjOhkTp1mli/MtDsZTw35g0o6VizV3FQr5U2lykolvBMqFQ9+DQr/vkIsX6ld6VL+SnHHoRdTEJepbVGCxM6qROv0fBIkoGmNRd0Pi5GTjyGO1gixdSQ/7mshGdCKd3eQWwONGK4URi48bKgw5g1dg3PmOlQl/qDA0FHkRcq1goTFpbwTDip4jy1h7qW4aAjmTE3JkhimgMDFpuQX8OzhGfCy03jJANs0suT/mYHuWRV0GEUXhG00lqjhTEmf0J+2dUSnjEmL8YbxzVsLOEZY/LHEp4xpa2izYUXXw46jMILuMtJLizhGVNA4kLlq2dxBweDDmV2hLyNyRKeMQUkadCdB4IOY9ZYCc+YEiRpaNrdh+w7jOuWyG1lYNfwjClFiV5Fn9sT9vM/vwLuVJwLS3jGFICkQ37mF4hVaY0pMaJQ/ej+sF+/LwxLeMZMT6Sxgb4FRXoParqErttlCPK2sVxYwjOhJdVVDNbZSCNFw67hGTNNIgwvbgg6CjMFQvgfR2IJz4STOHSuKw86CjNVVsIzprQ0/7YTt78/6DACYa20xkxD5JLV4a8fTUBKtMECsBKeMVMliQS9a+tQG562uGj4W2mzfqVEZImIPCEi+0Vkn4h80l/eICKPicgh/2e9v1xE5Ksi0iIiu0XkykL/EmbuiDQ20PbBK+lvtmxXlPI0xLuI3CsiZ0Rk7wTrp5VncvlWpYBPq+p64Brg4yKyHrgTeFxV1wCP+68B3g6s8aetwDdzCcQYAETQIq3Kmrw+texbwM2TrJ9Wnsma8FT1pKru9Od7gQPAIuBW4D5/s/uAd/nztwLfVs/TQJ2ILMglGGO6b1gTdAgzUn0shXvolaDDCE6eSniq+iugY5JNppVnplRvEJHlwBXAM8B8VT3przoFzPfnFwHHMt523F82dl9bRWSHiOxIMjSVMMwcFV2+lOHK4i7eOWlFU6mgwwjMFEp4TSPnvz9tneKhcsozY+XcaCEiVcAPgE+pao9kPGtTVVVkag3SqroN2AZQIw0hb9sxhRZdspi26xfhxoKOZPoiw1D2u6OUbButMpUBQNtUdXPhghlfTiU8EYnhJbv7VfUhf/HpkSKk//OMv/wEsCTj7Yv9ZcZMSCvKSBV5P+Py9jTp1tagwwjMyEN88nQNL5tp5ZlcWmkFuAc4oKp3ZazaDnzIn/8Q8KOM5R/0W1GuAbozqr7GzFmVvz0cdAjBm70HcU8rz+RSpX0D8AFgj4js8pd9FvgS8KCI3AEcBd7nr3sYuAVoAfqBP57Kb2FKjyQSdGxuCjqMGWk4MEi6oyvoMAInmp9sJiLfBd6Md63vOPB5IAagqnczzTyTNeGp6pNM3Of9xnG2V+DjuRzcGAARIVVR3I0VkYEkWkpDuY8nj6OlqOrtWdZPK8/YnRbGmLyxe2mNmYRTWYnTUB90GDPipMAZSJZu62yGsN9aZgnPBCq9aTVnLq0IOowZKe9Ik953MOgwwsFKeMaMT6JRXKe4r90BSOn2Mx4tf11OCsYSnglEdMFFnLl5BVrk30BJQ/ljL4S9YDN7Qv5BFPnXzRSj6OJFtL9pSVHfVTFKKY9/l2Gk43GYWcIzs8qprqb1xqWki/RhZGZy4oY741nCM7Nr1RJLdnOVPbXMmNE6N9QEHUJehb0KN9usW4oxGco600CEZJEPAzWi+denSZfwcFAXCPk/ABtH28yqxMPPUb+3O+gw8kaSluwyzeJoKdNiJTwz+w4d5aIznQAMXbyQnmXeRT03StE8uCfWrzQ+foR0R2fQoYSHAnkaPKBQLOGZWef29597bmvk5Cnq/cFknU2XMLCgkuGaCEO14a7yVp5Mkjp5KugwQseu4RmTjV8qcHftJ7ELKpqb0QWNpOrL6VhXFmxs46h9eZjYf1hn47GsH54x05BubYXWVuKLFhJfuBQANyakQpD7nBSUvdplDRXjUbUqrTHTlTrxGjXffQ3wbkUbWreQrlUJNBJcTFWvpUi/ZCMbT8RKeMbkQerkKSInTzH/lSXgOGg0Quv181Fnlhs6Ql6CCVzIPx5LeKaopI6efzJf4yuvAtD1h5tnpV9fZBgSj+8K+zkdqLCX8IqkE4AxF9JUCk2lqN/TRfXxwl9Tq20ZKOlnzmalQFpzmwJiJTxT9NzdL1L+UhkVtTV0vHUl6YTg5vmb7aQgdug1G9U4i7CX8CzhmTnBHRyEwUFq7z9DZO0qBlY20LM0f1/vhv1nSZ8+k33DUhfya5xWpTVzTvqlw5T9ck9+SxshH/YoLMJ+a5klPDMnuYODNN//Ag37B2bc+z/RrTiDyfwENpfl+hBuu5fWmPxz+/txntxF4/BlDCws5+y8qXfgExdqHtqJmxwuQIRziwASYINELqyEZ+a+Z/dQ+fgByjumXtRr3N2Hpqx0lytRzWkKiiU8UxLc3l4qXhucctKLdPSF/kJ8aBRBldYSnikZ8tQLVD+yn8oz1rmkMPT8/bTZpoBkTXgiskREnhCR/SKyT0Q+6S//goicEJFd/nRLxns+IyItInJQRG4q5C9gzFS4vb1U/mw35e3ZS3pOCiRlyXEq8tlKKyI3+zmkRUTuHGf9h0WkNSMH/Um2febSaJECPq2qO0WkGnheRB7z131FVf9hTBDrgduADcBC4OcislZV7ZtjQsEdHKT6Z3sZeP/GSberfWWI1JFXZymqOSJPpTcRiQBfB34POA48JyLbVXX/mE2/p6qfyHW/WUt4qnpSVXf6873AAWDRJG+5FXhAVYdU9RWgBdiSa0DGzAZ3YJCGg4OTlzbs0t3UqNdKm8uUgy1Ai6q+rKrDwAN4uWVGpnQNT0SWA1cAz/iLPiEiu0XkXhGp95ctAo5lvO044yRIEdkqIjtEZEeSoalHbsxMuGmcX++maoJ7cGP9SvTJ3bMc1ByQv0aLnPII8B4/B31fRJZk22nOCU9EqoAfAJ9S1R7gm8Aq4HLgJPDlXPcFoKrbVHWzqm6OYQ8qNQFw01Qc7sAZJ+fVHLaBAqZjCt1SmkYKPP60dRqH+3dguapuBB4D7sv2hpw6HotIDC/Z3a+qDwGo6umM9f8E/Nh/eQLIzLSL/WXGhE76pcPI6+eNOhMkDdF9r9hAAdOR+zW8NlXdPMn6rHlEVdszXv5f4O+zHTSXVloB7gEOqOpdGcsXZGz2bmCvP78duE1EEiKyAlgDPJvtOMYEpaJ1dEmu+flu0j09AUVTxBRwc5yyew5YIyIrRCSO1xC6PXODMTnonXjtC5PKpYT3BuADwB4R2eUv+yxwu4hcjvdrHgE+AqCq+0TkQWA/Xgvvx62F1oRZ+S/2wA2X0bs4SlmXIsdPZ3+TuYCQv7soVDUlIp8AHgEiwL1+bvlrYIeqbgf+XETeiZdnOoAPZ9tv1oSnqk/i3SY31sOTvOeLwBez7duYMHAHBylrHaR3cRWxs2nSbe3Z32TG5+bvOY2q+jBj8oyqfi5j/jPAZ6ayTxs8wBhf9bEU5Sd6c6xxmQuMVGlDzBKeMYBz6FXKdg/iDlkXqZkIcmCAXFjCMwZId3UHHcLcYAnPGFMa7EHcxphSMfLUshCzhGeMyRu7hmeMKR2W8IwxJUEJ/dPdLOEZY/LEGi2MMaXEEp4xpiQokA73rRaW8IwxeaKglvCMMaXCqrTGmJJgrbTGmJJiJTxjTMmwhGeMKQmqkA734OaW8Iwx+WMlPGNMybCEZ4wpDWqttMaYEqGg1vHYGFMy7NYyY0xJUM3rYxoLwRKeMSZ/rNHCGFMq1Ep4xpjSYAOAGmNKRREMHuBk20BEykTkWRF5QUT2icj/9JevEJFnRKRFRL4nInF/ecJ/3eKvX17g38EYEwIKaDqd0xSUrAkPGAJuUNVNwOXAzSJyDfB3wFdUdTXQCdzhb38H0Okv/4q/nTFmrlN/ANBcphyIyM0ictAvPN05zvopF66yJjz19PkvY/6kwA3A9/3l9wHv8udv9V/jr79RRCTbcYwxxU9dzWnKRkQiwNeBtwPrgdtFZP2YzaZcuMqlhIeIRERkF3AGeAw4DHSpasrf5DiwyJ9fBBwD8Nd3A43j7HOriOwQkR1JhnIJwxgTdvkr4W0BWlT1ZVUdBh7AK0xlmnLhKqdGC1VNA5eLSB3wQ+DiXN6XZZ/bgG0AItL6c/3+WaBtpvvNoyYsnslYPJMrxniWzeQAvXQ+8nP9flOOm5eJyI6M19v8nDDiXMHJdxy4esw+RhWuRGSkcDXh7zmlVlpV7RKRJ4DXA3UiEvVLcYuBE/5mJ4AlwHERiQK1QHuW/TaLyA5V3TyVeArJ4pmcxTO5UoxHVW8u5P7zIZdW2ma/ZIeIlAO/BxwAngDe62/2IeBH/vx2/zX++l+ohrxzjjEmbEYKTiMyC1UXbJNr4SqXa3gLgCdEZDfwHPCYqv4Y+B/AX4pIC14x8h5/+3uARn/5XwIXtK4YY0wWzwFr/O5vceA2vMJUpikXrrJWaVV1N3DFOMtfxruwOHb5IPAH2fY7jm3ZN5lVFs/kLJ7JWTwz4F+T+wTwCBAB7lXVfSLy18AOVd2OV7j6F79w1YGXFCclVts0xpSKnLqlGGPMXGAJzxhTMgJPeNluH5mlGI6IyB4R2TXSN0hEGkTkMRE55P+sL+Dx7xWRMyKyN2PZuMcXz1f9z2u3iFw5izF9QURO+J/TLhG5JWPdZ/yYDorITXmOZYmIPCEi+/37uT/pLw/kM5oknkA+H3//ds97LlQ1sAnvYuRhYCUQB14A1gcQxxGgacyyvwfu9OfvBP6ugMd/I3AlsDfb8YFbgJ8CAlwDPDOLMX0B+G/jbLve/9slgBX+3zSSx1gWAFf689XAS/4xA/mMJoknkM/HP4YAVf58DHjG/90fBG7zl98NfNSf/xhwtz9/G/C9Qn2/wzQFXcLL5faRoGTetpJ5r3Deqeqv8FqZcjn+rcC31fM0XgfwBbMU00RuBR5Q1SFVfQVoYZwW/BnEclJVd/rzvXj9QBcR0Gc0STwTKejn48ehave8ZxV0whvv9pHJvjiFosCjIvK8iGz1l81X1ZP+/Clg/izHNNHxg/7MPuFXE+/NqObPWkx+1esKvBJM4J/RmHggwM9HCnDP+1wTdMILi+tU9Uq8kRk+LiJvzFypXrk/sP47QR8/wzeBVXjDhJ0EvjybBxeRKuAHwKdUtSdzXRCf0TjxBPr5qGpaVS/HuythC3m4532uCTrh5XL7SMGp6gn/5xm8wRG2AKdHqkH+zzOzHNZExw/sM1PV0/5J5QL/xPlqWcFjEpEYXnK5X1Uf8hcH9hmNF0+Qn08mVe3Cu/Xz3D3v4xx3yrdlzQVBJ7xcbh8pKBGpFJHqkXngbcBeRt+2knmv8GyZ6PjbgQ/6LZHXAN0Z1bqCGnMd7N14n9NITLf5LX8rgDXAs3k8ruD1qj+gqndlrArkM5oonqA+H//Yds97LoJuNcFrUXsJ73rDXwVw/JV4LWgvAPtGYsC7nvE4cAj4OdBQwBi+i1cFSuJdZ7ljouPjtcZ93f+89gCbZzGmf/GPuRvvhFmQsf1f+TEdBN6e51iuw6uu7gZ2+dMtQX1Gk8QTyOfj738j8Dv/2HuBz2V8v5/Fayj5VyDhLy/zX7f461cW6vsdpsluLTPGlIygq7TGGDNrLOEZY0qGJTxjTMmwhGeMKRmW8IwxJcMSnjGmZFjCM8aUjP8PkO6W4f1Oca0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# verification\n",
"\n",
"tmp = np.zeros(llung.shape)\n",
"tmp[llung==1]=1\n",
"tmp[nonfiss==1]=2\n",
"tmp[fiss==1]=3\n",
"\n",
"ind = 200\n",
"plt.imshow(tmp[:,:,ind].squeeze(),vmin=0,vmax=3)\n",
"plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#\n",
"# nonfiss==1\n",
"# fiss==1\n"
]
},
{
"cell_type": "code",
"execution_count": 120,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"89 240 205\n",
"90 239 203\n",
"90 239 204\n",
"90 239 205\n",
"90 240 203\n",
"90 240 204\n",
"90 240 205\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOaUlEQVR4nO3db4hc53XH8d/ZXe2f0X+tVrIqCUe2ZNlGLk6xQt2Ktq7SojihCqYUG1KcNpBXadNSCAp94beBlpJCQ4tJXBnq2hQ3JqamioTbWARaE1kRjWU5lVFVWWv9Wwu8Wku78u6evtipO1lL1fPMvXPvnPH3A2ZnZo+fOTN7OHt19555zN0FAIinr+4EAADtoYEDQFA0cAAIigYOAEHRwAEgKBo4AAQ1UOWTDQ0NeaPRSIodGEhPLedSyGvXriXHzs7OdiSH+fn50tfslNRcc5lZcmxfX9pxxvz8vObn59MXLpGZ1f/DQk9z94/UdqUNvNFo6KGHHkqKXbduXfK609PTybHHjx9Pjp2YmEiOzWn2U1NTpa+ZI6cp5/zCy/mFMzQ0lBw7PDycFHflypXkNYFeUOgUipntMbOfmtlbZravrKSAulHbiKDtBm5m/ZK+Jekzku6V9JiZ3VtWYkBdqG1EUeQI/FOS3nL3U+5+XdJzkvaWkxZQK2obIRRp4Bslvd1y/2zzsZ9hZl82syNmdmRmZqbA0wGVya7tyjIDWnT8MkJ3f9LdH3D3B3L+cAV0u9barjsXfDwVaeDjkja33N/UfAyIjtpGCEUa+I8kbTOzLWY2KOlRSS+WkxZQK2obIbR9Hbi7z5rZVyR9X1K/pKfcPf0ia6BLUduIwqqc9hsdHfWHH344KXZsbCx53TvuuCM5NmeI5cCBA8mxJ06cSI5NHTjJGVDK+Tn29/d3ZN25ubnk2MHBwdJzmJqa0tzcHJOY6Ek3msTks1AAICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABFXpnpirV6/WI488khR7/vz55HVzRri3bNmSHPvggw8mx46Pl/9hdakbQEudG4/P+eiB1L0+Jen69evJsal7YuZslAz0Ao7AASCoIntibjazfzWzN8zsuJl9tczEgLpQ24iiyCmUWUl/4u5HzWy5pNfM7JC7v1FSbkBdqG2E0PYRuLufc/ejzdtXJJ3QDfYNBKKhthFFKefAzewTkj4p6dUy1gO6BbWNbla4gZvZMkn/KOmP3H3yBt//cOfuycmPfBvoWjm1XX12QMEGbmZLtFDgz7j7d28U07pz94oVK4o8HVCZ3NquNjtgQZGrUEzSdySdcPe/KC8loF7UNqIocgT+y5J+V9Kvm9mx5n9pG14C3Y3aRghFdqX/oSRG39BzqG1EUekofaPR0M6dO5NiJyYmktfN2RH+3LlzybH33XdfcuypU6eSY1966aWkuCVLliSvmTNKnzPyPjs7mxyb85EGQ0NDybHXrl1Lisv5iACgFzBKDwBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoSkfpcwwMpKc2NjaWHJszcp6TQ84O9seOHUuKu3jxYvKaOfr6OvN7e2RkJDk2Z5Q+9eeQs9M90As4AgeAoMrYkaffzH5sZv9URkJAt6C20e3KOAL/qhY2fQV6DbWNrlZ0S7VNkj4r6dvlpAN0B2obERQ9Av+mpK9Jmi+eCtBVvilqG12uyJ6Yn5N00d1fu0Xchzt3X758ud2nAyrTTm1XlBrwM4ruiflbZnZa0nNa2D/w7xYHte7cvWbNmgJPB1Qmu7arThCQCjRwd/+6u29y909IelTSv7j7F0rLDKgJtY0ouA4cAIIqZRLT3X8g6QdlrAV0E2ob3azSUXozS95pfdmyZcnrbtq0KTk2Z92JiYnk2JUrVybH7tq1Kynu8OHDyWvm7DSfs9t9o9FIjjWz5Ngco6OjSXFXrlzpyPMD3YpTKAAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCqnSUfmBgIHksenZ2Nnldd0+Ove2225JjT548mRx75Ej6R0KvW7cuKW7z5s3Jax49ejQ5dnp6Ojk2Z+x+1apVybFjY2PJsVu3bk2Ke+edd5LXBHoBR+AAEFTRPTFXmdnzZvammZ0wswfLSgyoE7WNCIqeQvlLSQfc/bfNbFBS+kfXAd2N2kbXa7uBm9lKSb8i6YuS5O7XJV0vJy2gPtQ2oihyCmWLpEuS/tbMfmxm3zazpYuDWjd+zfl8baBG2bVdfYpAsQY+IOkXJP21u39S0vuS9i0Oat34de3atQWeDqhMdm1XnSAgFWvgZyWddfdXm/ef10LRA9FR2wihyK705yW9bWbbmw/tlvRGKVkBNaK2EUXRq1D+QNIzzb/Sn5L0e8VTAroCtY2uV6iBu/sxScnn/+bm5jQ5OVnkKW+6bqq+vvR/dCxfvjw59vTp08mxr7zySlLcyMhI8pp33nlncmzOH5Nzpiu3b99+66CmHTt2JMfeddddSXEHDx5MXvNWcmsbeXKmpzu1WXYvYBITAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCooEDQFCVbmrc39+vRiNtY5OrV68mr5szcj48PJwcmzPuu3HjxuTYgYG0tz1n5P2ee+5Jjt29e3dybM6I/t13350ce/vttyfHrlixIilu5cqVyWuiXozHl4MjcAAIigYOAEEV3ZX+j83suJm9bmbPmln6+Qmgi1HbiKDtBm5mGyX9oaQH3H2HpH5Jj5aVGFAXahtRFD2FMiBpxMwGJDUkvVM8JaArUNvoekW2VBuX9OeSzkg6J+k9d//IJ+q37tx96dKl9jMFKtJObVedIyAVO4WyWtJeSVsk/ZykpWb2hcVxrTt3j42NtZ8pUJF2arvqHAGp2CmUT0v6L3e/5O4fSPqupF8qJy2gVtQ2QijSwM9I+kUza9jCVfm7JZ0oJy2gVtQ2QihyDvxVSc9LOirpJ821niwpL6A21DaiKLor/ROSnkiNn5+f1/T0dFLs+++/n5xH6pqSNDMzkxw7OzubHJszSp86cj4+Pp68Zn9/f3Ls1q1bk2N37tyZHLthw4bk2Jz39sKFC0lxH3zwQfKat5Jb20AdmMQEgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEFSlu9LPzc1pamoqKXZwcDB53SVLliTH5oxbL1u2LDl2165dybFr1qxJijtz5kzymvPz88mxIyMjybE571fOz6yvL/3YgR3MgRvjCBwAgrplAzezp8zsopm93vLYGjM7ZGYnm19XdzZNoHzUNqJLOQLfL2nPosf2SXrZ3bdJerl5H4hmv6htBHbLBu7uhyVdXvTwXklPN28/Lenz5aYFdB61jejaPQe+3t3PNW+fl7S+pHyAulHbCKPwHzHd3SX5zb7funP35cuLD3aA7pVT2xWmBXyo3QZ+wcw2SFLz68WbBbbu3J16+RxQo7Zqu7LsgBbtNvAXJT3evP24pO+Vkw5QO2obYaRcRvispH+TtN3MzprZlyR9Q9JvmNlJSZ9u3gdCobYR3S0nMd39sZt8a3fJuQCVorYRXaWj9J2SM5adM3a/dOnS5Njh4eHk2Pvvvz8pbtu2bclrTk9PJ8cODKT/2HN2u5+YmEiOzRnnT43NqQOgF1DxABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKqdJS+r68veeQ8Z+Q9J3Zqaio59r333kuOHRoaSo5NtWLFiuTYnNH0Tr23V69eTY7Neb9Sc8gZ+wd6AUfgABBUu7vS/5mZvWlm/2FmL5jZqo5mCXQAtY3o2t2V/pCkHe7+85L+U9LXS84LqMJ+UdsIrK1d6d39oLvPNu/+u6RNHcgN6ChqG9GVcQ789yX9882+2brx67vvvlvC0wGVSa7tCnMCPlSogZvZn0qalfTMzWJaN34dHR0t8nRAZXJru7rMgP/T9mWEZvZFSZ+TtNvdvbSMgJpR24iirQZuZnskfU3Sr7p7+sW/QJejthFJu7vS/5Wk5ZIOmdkxM/ubDucJlI7aRnTt7kr/nQ7kAlSK2kZ0lY7Sz83NaXJyMik2Z4zczJJjc0a4c05/5oyRz8zMJMU1Go3kNXPGyHPer5yd3nPG7gcG0ksv9WeW87qAXsAoPQAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgrMpPyzSzS5L+e9HDayVNVJZEtXr1tXXr67rd3cfqeOKPWW336uuSuve13bC2K23gN2JmR3r1A/F79bX16usqW6++T736uqR4r41TKAAQFA0cAILqhgb+ZN0JdFCvvrZefV1l69X3qVdflxTstdV+DhwA0J5uOAIHALSh1gZuZnvM7Kdm9paZ7aszlzKZ2Wkz+0lzT8UjdedThJk9ZWYXzez1lsfWmNkhMzvZ/Lq6zhy7EbXd/Xqhtmtr4GbWL+lbkj4j6V5Jj5nZvXXl0wEPufv9kS5Juon9kvYsemyfpJfdfZukl5v30URth7FfwWu7ziPwT0l6y91Puft1Sc9J2ltjPrgBdz8s6fKih/dKerp5+2lJn68ypwCo7QB6obbrbOAbJb3dcv9s87Fe4JIOmtlrZvblupPpgPXufq55+7yk9XUm04Wo7bhC1Xalu9J/jOxy93EzWyfpkJm92fxt33Pc3c2MS5k+PqjtLlLnEfi4pM0t9zc1HwvP3cebXy9KekEL/6TuJRfMbIMkNb9erDmfbkNtxxWqtuts4D+StM3MtpjZoKRHJb1YYz6lMLOlZrb8f29L+k1Jr////1c4L0p6vHn7cUnfqzGXbkRtxxWqtms7heLus2b2FUnfl9Qv6Sl3P15XPiVaL+kFM5MW3t+/d/cD9abUPjN7VtKvSVprZmclPSHpG5L+wcy+pIVP4Pud+jLsPtR2DL1Q20xiAkBQTGICQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgvofDtc8A1GoWDcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOaklEQVR4nO3de4yc5XXH8d/xXn1j18Z0fY8tMBEQRbiygtOiXsCN7DaKI1FVICWQi5S/cmlVKXKUP/g3UlCVSolaoZgaqRRU0SCjSiVYtFFUqSBfYiWYJYAhiXdjYhvLGLw37+7pHzshk40dnmfed953zvj7kaKdGZ88c97Zw/H4nfeZY+4uAEA8S+pOAADQGho4AARFAweAoGjgABAUDRwAgqKBA0BQvVU+WX9/vw8ODibFLlmS/nfLzMxMcuzly5eTY3t6epJjcy7HnJ+fr21NKe+4cn4POa9tzrGl5jA/Py93t+SFS2RmXI+LtrpSbVfawAcHB7Vz587k2FSnTp1Kjh0fH0+OHRoaSo6dnZ1Njp2YmEiKm56eTl4z5y+x6667Ljl2YGAgOfb06dPJsXNzc8mxqbUwNTWVvCbQDQqdQjGz3Wb2UzN7zcz2lZUUUDdqGxG03MDNrEfSdyTtkXSrpPvM7NayEgPqQm0jiiLvwD8i6TV3f93dZyQ9IWlvOWkBtaK2EUKRBr5BUvPJ57HGY7/FzL5gZkfM7EjOh1xAjbJru7LMgCZtv4zQ3R929x3uvqOvr6/dTwdUprm2684F16YiDXxc0qam+xsbjwHRUdsIoUgDPyxpm5ltNbN+SfdKerqctIBaUdsIoeXrwN191sy+KOn7knokPeLuJ0rLDKgJtY0orMqBDitXrvQdO9JOFy5dujR53TVr1iTHXrhwITn2/PnzbVk3dcNJzuacycnJ5Nic37lZ+sbGnI1H7dh0Mzs7q/n5eXZioitdaScm34UCAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABIKhKZ2IODw9r796078U/fvx48ro5A31vvPHG5NhbbrklOTZnUPDzzz+fFHfy5MnkNXO20rdrdmTOAOScmZg5v1/gWsI7cAAIqshMzE1m9j9m9pKZnTCzr5SZGFAXahtRFDmFMivp7939mJmtlHTUzA65+0sl5QbUhdpGCC2/A3f30+5+rHH7HUmjusLcQCAaahtRlHIO3My2SNou6YUy1gM6BbWNTla4gZvZCkn/Ielv3f3iFf78vcndly5dKvp0QGVyarv67ICCDdzM+rRQ4I+5+/euFNM8uXv58uVFng6oTG5tV5sdsKDIVSgmab+kUXf/h/JSAupFbSOKIu/A/1jSpyXdZWbHG//7y5LyAupEbSOEIlPp/1dSLQNkgXaithFFpVvph4aGtGfPnqTYkZGR5HVHR0eTY3O2nK9YsSI59rbbbkuOTd3KPj4+nrxmzuuV8xrkbGPP2R6fIzXfnOMCugFb6QEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAEVelW+r6+Pq1duzYp9s4770xed/369cmxR48eTY49ciT9a57HxsaSY2+66aakuC1btiSveeLEieTYgYGB5Ngcb7/9dnKsuyfH5ky7B64l/JcBAEGVMZGnx8x+ZGb/WUZCQKegttHpyngH/hUtDH0Fug21jY5WdKTaRkl/Jem75aQDdAZqGxEUfQf+LUlflZT+pdFADN8StY0OV2Qm5sclnXH333tZR/Pk7rfeeqvVpwMq00ptV5Qa8FuKzsT8hJn9TNITWpgf+K+Lg5ond19//fUFng6oTHZtV50gIBVo4O7+NXff6O5bJN0r6b/d/VOlZQbUhNpGFFwHDgBBlbIT091/IOkHZawFdBJqG52s0q30c3NzunjxYunrbt68OTm2p6cnOTZnInvOBPlUd911V3JszgfEOVvec7bdv/POO8mxOVvp+/r6kuKYSo9rDadQACAoGjgABEUDB4CgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARV6VZ6STKzpLicLe/9/f3JsSMjI8mxO3fuTI4dHU2fvPXKK68kxQ0NDSWvuX379uTYw4cPJ8e+++67ybGpW94lqbe3/NJLrS2gW/AOHACCKjoTc9jMnjSzl81s1Mw+WlZiQJ2obURQ9N+x/yjpGXf/azPrl7SshJyATkBto+O13MDNbEjSn0j6jCS5+4ykmXLSAupDbSOKIqdQtko6K+lfzOxHZvZdM1u+OKh58Ov58+cLPB1Qmezarj5FoFgD75X0h5L+yd23S7okad/ioObBr6tXry7wdEBlsmu76gQBqVgDH5M05u4vNO4/qYWiB6KjthFCkan0b0o6ZWYfbDx0t6SXSskKqBG1jSiKXoXyJUmPNT6lf13SZ4unBHQEahsdr1ADd/fjkpLP/y1ZskTLlqVdjZWzU2/JkvR/SOSsu3z573xudVXDw8PJsakDiFN3bOY+/80335wce+zYseTYnCHQOb+z1Ngyd2Lm1jby5Ay1bpdu2LnLTkwACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQVU61Li3t1erVq1Kip2dnU1ed2pqKjk2Z1jy4OBgcmzOsOTU2IMHDyavmfNd6+vWrUuO3bRpU3Ls2NhYcuzly5eTY1esWJEUl7M9H8jZzt+p2+6peAAIigYOAEEVnUr/d2Z2wsxeNLPHzSz9nAPQwahtRNByAzezDZK+LGmHu39IUo+ke8tKDKgLtY0oip5C6ZW01Mx6JS2T9MviKQEdgdpGxysyUm1c0kOSfiHptKS33f3ZxXHNk7vPnj3beqZARVqp7apzBKRip1BWSdoraauk9ZKWm9mnFsc1T+6+4YYbWs8UqEgrtV11joBU7BTKLklvuPtZd78s6XuS/qictIBaUdsIoUgD/4WknWa2zBaucr9b0mg5aQG1orYRQpFz4C9IelLSMUk/aaz1cEl5AbWhthFF0an0D0p6MCNe09PTSbE529hzpqFPTk4mx+Zszc7Zop+6lT3nNZiYmEiOPXnyZHLsxo0bk2PvuOOO5Ni1a9cmx27evDkp7qGHHkpe8/3k1jbydMLW9Jyt9J267Z6dmAAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCqnQqvZmpv78/KTZn62rONva5ubnk2Jyt9KnHJUnbtm1LituxI/1bSi9cuJAcmzPB/ty5c8mxW7duTY69/fbbk2N37dqVFLd///7kNYGcLe85/ahKvAMHgKDet4Gb2SNmdsbMXmx6bLWZHTKzVxs/V7U3TaB81DaiS3kHfkDS7kWP7ZP0nLtvk/Rc4z4QzQFR2wjsfRu4u/9Q0uKTpnslPdq4/aikT5abFtB+1Daia/Uc+Ii7n27cflPSSEn5AHWjthFG4Q8xfeHj2at+RMtUekSVU9sVpgW8p9UG/iszWydJjZ9nrhbIVHoE01JtV5Yd0KTVBv60pAcatx+QdLCcdIDaUdsII+Uywscl/Z+kD5rZmJl9XtI3JP2Fmb0qaVfjPhAKtY3o3ncnprvfd5U/urvkXIBKUduIrtKt9PPz85qamkqKzdmaPjMzkxybs+0+Z9r99PR0cuz69euT4u6///7kNTds2JAc+8YbbyTHLl26tC2xy5YtK33dnK8+AHJUOWk+BxUPAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABIKjKp9L39fUlxU5OTiavOzs7mxw7ODiYHNuubfepOeRMeb/nnnuSYycmJpJjU7/6QMr7PQwMDCTHXrx4MSlubm4ueU2gG/AOHACCanUq/TfN7GUz+7GZPWVmw23NEmgDahvRtTqV/pCkD7n7hyW9IulrJecFVOGAqG0E1tJUend/1t1/fcLzeUkb25Ab0FbUNqIr4xz45yT919X+sHnw67lz50p4OqAyybVdYU7Aewo1cDP7uqRZSY9dLaZ58OuaNWuKPB1Qmdzari4z4DdavozQzD4j6eOS7nZ3Ly0joGbUNqJoqYGb2W5JX5X0p+6eflEx0OGobUTS6lT6b0taKemQmR03s39uc55A6ahtRNfqVPr9bcgFqBS1jegq3Uo/NzeXvC06Z8J4zhbqS5cuJcfmyMkhdYJ9znbz4eHhtsTmyJncnfqVCpI0MzOTFNfbW2k5A7VjKz0ABEUDB4CgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoKzKb8s0s7OSfr7o4TWSunXSQ7ceW6ce1wfc/YY6nvgaq+1uPS6pc4/tirVdaQO/EjM70q1fiN+tx9atx1W2bn2duvW4pHjHxikUAAiKBg4AQXVCA3+47gTaqFuPrVuPq2zd+jp163FJwY6t9nPgAIDWdMI7cABAC2pt4Ga228x+amavmdm+OnMpk5n9zMx+0pipeKTufIows0fM7IyZvdj02GozO2RmrzZ+rqozx05EbXe+bqjt2hq4mfVI+o6kPZJulXSfmd1aVz5t8OfufnukS5Ku4oCk3Yse2yfpOXffJum5xn00UNthHFDw2q7zHfhHJL3m7q+7+4ykJyTtrTEfXIG7/1DS+UUP75X0aOP2o5I+WWVOAVDbAXRDbdfZwDdIOtV0f6zxWDdwSc+a2VEz+0LdybTBiLufbtx+U9JIncl0IGo7rlC1zRjv9rjT3cfN7A8kHTKzlxt/23cdd3cz41Kmawe13UHqfAc+LmlT0/2NjcfCc/fxxs8zkp7Swj+pu8mvzGydJDV+nqk5n05DbccVqrbrbOCHJW0zs61m1i/pXklP15hPKcxsuZmt/PVtSR+T9OLv/3+F87SkBxq3H5B0sMZcOhG1HVeo2q7tFIq7z5rZFyV9X1KPpEfc/URd+ZRoRNJTZiYtvL7/5u7P1JtS68zscUl/JmmNmY1JelDSNyT9u5l9XgvfwPc39WXYeajtGLqhttmJCQBBsRMTAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAENT/A9H7bAxh2QhpAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOg0lEQVR4nO3db4xc5XXH8d/Z9a7X68VrzMLiv7UBE2GiCldWoAXRFpzitFEdiaoCKYW0kfwqbVpVihz1BW8jtapSqVErlFAjlYIqGmRUqQkWbRRVKshrBzWYJYAhjXdr4jUrY8P+t09f7CSdbOzwPHPv3Dtn/P1IaGfHh+ee2T06vr5znznm7gIAxNNTdwIAgNbQwAEgKBo4AARFAweAoGjgABAUDRwAglpV5cH6+/t9cHAwKXZhYSF53cXFxeTY3t7e5NicWywvXbpU+ro5a+a8rp6e9L+3c362OT+vnBzMLCnu0qVLunTpUlpwycyM+3HRVu7+c7VdaQMfHBzUvffemxR76tSp5HUnJyeTY4eHh5Njl5aWkmNnZmaSY+fn55Picv4SW7duXXLs6tWrk2NPnz6dHHvx4sXk2IGBgeTY/v7+pLgLFy4krwl0g0KXUMxsn5n9wMzeMrODZSUF1I3aRgQtN3Az65X0NUmfkrRL0sNmtqusxIC6UNuIosgZ+CckveXub7v7gqRnJO0vJy2gVtQ2QijSwDdLar5QPdF47meY2QEzGzOzsZxrukCNsmu7ssyAJm2/jdDdH3f3Pe6+J/XNKCCC5tquOxdcnYo08ElJW5u+39J4DoiO2kYIRRr4UUk7zWyHmfVLekjS8+WkBdSK2kYILd8H7u5LZvYFSd+W1CvpCXc/UVpmQE2obURhVQ50GB4e9rvvvjspdmRkJHndc+fOJcdOT0+3Zd25ubnk2NQ3c2dnZ5PXzPk9pu5slNI3HUl5P4McfX19yce/ePEiOzHRlS63E5PPQgGAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQVKUzMfv6+nTDDTckxeYM9L355puTY2+77bbk2JxBwS+99FJy7MmTJ5PicrbSt2sbe87w4ZyZmHUPgQa6AWfgABBUkZmYW83sP8zsNTM7YWZfLDMxoC7UNqIocgllSdKfu/txM7tG0jEzO+Lur5WUG1AXahshtHwG7u6n3f144/EFSeO6zNxAIBpqG1GUcg3czLZL2i3p5TLWAzoFtY1OVvguFDMbkvQvkv7U3c9f5s8PSDogSWvXri16OKAyObUN1KHQGbiZ9Wm5wJ9y929eLqZ5cvfAwECRwwGVya3tarMDlhW5C8UkfUPSuLv/dXkpAfWithFFkTPwuyX9gaT7zOyVxn+/XVJeQJ2obYRQZCr9f0qqZYAs0E7UNqKodCv9unXr9MADDyTFjo+PJ6+bs+V8aGgoOfb2229Pjs3Zyj45OZkUNzo6mrxmzs8gZ8t5zvb4HDn5Li4uJsUtX/kArh5spQeAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQVKVb6YeGhnTPPfckxW7atCl53WPHjiXHjo2NJcdOTEwkx95yyy3Jsdu3b0+KO3HiRPKaq1evTo7N8f777yfHpk6Pl/Km3ad+DPHS0lLymkA34AwcAIIq3MDNrNfMvmdm/1pGQkCnoLbR6co4A/+iloe+At2G2kZHKzpSbYuk35H09XLSAToDtY0Iip6Bf1XSlySlf8A0EMNXRW2jwxWZiflpSWfc/RfeAmJmB8xszMzGpqenWz0cUJlWarui1ICfUXQm5u+a2Q8lPaPl+YH/uDKoeXL3hg0bChwOqEx2bVedICAVaODu/mV33+Lu2yU9JOnf3f2zpWUG1ITaRhTcBw4AQZWyE9PdvyPpO2WsBXQSahudrNKt9Dm2bduWHNvb25scmzORPXV6fK777rsvKe69995LXjNny3vOtvsLFy4kx+Zspe/r60uOTZ02PzMzk7wm0A24hAIAQdHAASAoGjgABEUDB4CgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgqEq30puZVq1KO2R/f3/yuqOjo8mxd911V3Ls+Hj6NK033ngjOXZ4eDgpbvfu3clrHj16NDn2gw8+SI7N2fKe+rvNNT8/35Z1geg4AweAoIrOxFxvZs+a2etmNm5mv1pWYkCdqG1EUPTfvH8j6Vvu/ntm1i9psIScgE5AbaPjtdzAzWxY0r2SPidJ7r4gaaGctID6UNuIosgllB2SpiT9g5l9z8y+bmZrVwY1D37N+XxroEbZtV19ikCxBr5K0q9I+jt33y3pQ0kHVwY1D3697rrrChwOqEx2bVedICAVa+ATkibc/eXG989queiB6KhthFBkKv27kk6Z2ccaT90v6bVSsgJqRG0jiqJ3ofyxpKca79K/LekPi6cEdARqGx2vUAN391ckJV//6+np0dq1P/de0BVjU+XsAEw9viStX78+OTbnDdrUXZs5x7/11luTY48fP54cmzMEOud3lhM7MDCQFHf+/PnkNT9Kbm1jWc5g63ZIHYDdLdiJCQBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoSoca9/T0aM2aNUmxc3Nzyev29vYmx6Zuy5byhiXnxB4+fDgpbnp6OnnNjRs3Jsdu3bo1OXZiYiI5dnFxMTl2aGgoOfamm25KipuamkpeE90pZyt/N2y75wwcAIKigQNAUEWn0v+ZmZ0ws1fN7GkzS78+AXQwahsRtNzAzWyzpD+RtMfdPy6pV9JDZSUG1IXaRhRFL6GskrTGzFZJGpT0v8VTAjoCtY2OV2Sk2qSkv5L0I0mnJb3v7i+sjGue3H327NnWMwUq0kptV50jIBW7hHKtpP2SdkjaJGmtmX12ZVzz5O6RkZHWMwUq0kptV50jIBW7hLJX0jvuPuXui5K+KenXykkLqBW1jRCKNPAfSbrLzAZt+Y74+yWNl5MWUCtqGyEUuQb+sqRnJR2X9P3GWo+XlBdQG2obURSdSv+YpMdS481MfX19SbE509BnZ2eTY3Omoeds0c/Zyp66nX9mZiZ5zZMnTybHbtmyJTn2zjvvTI698cYbk2O3bduWHLtr166kuPHx8k6Sc2sby+renp6zlb4btt2zExMAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUJVOpZfSt6/mbGO/ePFicmzOVvr+/v7k2J07dybH7tmT9umj586dS14zZ4J9zuey79ixIzn2jjvuSI7du3dvcmzqBPucSffoTjlb3nO20ncqzsABIKiPbOBm9oSZnTGzV5ue22BmR8zszcbXa9ubJlA+ahvRpZyBH5K0b8VzByW96O47Jb3Y+B6I5pCobQT2kQ3c3b8raeUF1v2Snmw8flLSZ8pNC2g/ahvRtXoNfNTdTzcevytptKR8gLpR2wij8JuYvvxW7hXfzm2e3D01NVX0cEBlcmq7wrSAn2q1gf/YzDZKUuPrmSsFNk/uvv7661s8HFCZlmq7suyAJq028OclPdp4/Kikw+WkA9SO2kYYKbcRPi3pvyR9zMwmzOzzkr4i6ZNm9qakvY3vgVCobUT3kTsx3f3hK/zR/SXnAlSK2kZ0lW6ld/fkbe8LCwvJ6+Zsu8+Zdj8/P58cu2nTpuTYRx55JClu8+bNyWu+8847ybFr1qxpS+zg4GBb1u2GLc/oPJ06aT4HW+kBICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABFX5VvrZ2dmk2KWlpeR1BwYGkmPbte0+J4fUSe8PPvhg8pozMzPJsXNzc8mxOb+H1atXJ8eeP38+OTZ1i37O7wvoBpyBA0BQrU6l/0sze93M/tvMnjOz9W3NEmgDahvRtTqV/oikj7v7L0t6Q9KXS84LqMIhUdsIrKWp9O7+grv/5OLoS5K2tCE3oK2obURXxjXwP5L0b1f6w+bBr2fPni3hcEBlkmu7wpyAnyrUwM3sLyQtSXrqSjHNg19HRkaKHA6oTG5tV5cZ8P9avo3QzD4n6dOS7ndGpqCLUNuIoqUGbmb7JH1J0q+7e/oNyECHo7YRSatT6f9W0jWSjpjZK2b2923OEygdtY3oWp1K/4025AJUitpGdJVvpU/dmp06vV6SPvzww1ZTKi2HnAn2qVvO169fn7xmTmyOnMndfX19ybELCwulr5vzMQlAN2ArPQAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgrMpPyzSzKUn/s+LpEUndOumhW19bp76uX3L36+s48FVW2936uqTOfW2Xre1KG/jlmNlYt34gfre+tm59XWXr1p9Tt74uKd5r4xIKAARFAweAoDqhgT9edwJt1K2vrVtfV9m69efUra9LCvbaar8GDgBoTSecgQMAWlBrAzezfWb2AzN7y8wO1plLmczsh2b2/cZMxbG68ynCzJ4wszNm9mrTcxvM7IiZvdn4em2dOXYiarvzdUNt19bAzaxX0tckfUrSLkkPm9muuvJpg9909zsi3ZJ0BYck7Vvx3EFJL7r7TkkvNr5HA7UdxiEFr+06z8A/Iektd3/b3RckPSNpf4354DLc/buSplc8vV/Sk43HT0r6TJU5BUBtB9ANtV1nA98s6VTT9xON57qBS3rBzI6Z2YG6k2mDUXc/3Xj8rqTROpPpQNR2XKFqu9Kp9FeRe9x90sxukHTEzF5v/G3fddzdzYxbma4e1HYHqfMMfFLS1qbvtzSeC8/dJxtfz0h6Tsv/pO4mPzazjZLU+Hqm5nw6DbUdV6jarrOBH5W008x2mFm/pIckPV9jPqUws7Vmds1PHkv6LUmv/uL/K5znJT3aePyopMM15tKJqO24QtV2bZdQ3H3JzL4g6duSeiU94e4n6sqnRKOSnjMzafnn+0/u/q16U2qdmT0t6TckjZjZhKTHJH1F0j+b2ee1/Al8v19fhp2H2o6hG2qbnZgAEBQ7MQEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB/R/oLmw4dfuOcgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOdklEQVR4nO3dbYxU93XH8d+PXXZhARsH3DWPBRuMhFFlKhS7D+oTtCJ1FPKiqmwptdNG4lXatKoUEfWF30ZqVaVSo1ZW4mKprq3KjYVdqYmR2yiqFFs8xGqMDbax0xiMAxhBcXnYXfb0xU7SDYbw/8+9c+/8h+9HsnZmOL5z7nB09nLnnvt3RAgAUJ45bScAAOgODRwACkUDB4BC0cABoFA0cAAoFA0cAAo13OibDQ/HyMhIUuzk5GTydoeGhpJjcy6bnJ6ebnW7Ofs1Z0767+KczzZnv3JysF37dqempjQ9PZ2+4RrZ5npc9FREfKS2G23gIyMj2rBhQ1Lse++9l7zdW265JTl2amoqOfbChQvJsZcvX06OnZiYSIrL2a/R0dHk2BMnTiTHXrlyJTl23rx5ybGpv8glaf78+Ulxp0+fTt4mMAgqnUKxvd32Edtv2d5VV1JA26htlKDrBm57SNJXJX1C0kZJD9neWFdiQFuobZSiyhH4xyW9FRFvR8SEpKcl7agnLaBV1DaKUKWBr5D07qznxzqv/RTbO23vt70/5/wz0KLs2m4sM2CWnl9GGBGPRcSWiNgyPNzod6ZAT82u7bZzwc2pSgM/LmnVrOcrO68BpaO2UYQqDXyfpPW219oekfSgpOfqSQtoFbWNInR9TiMipmx/XtK3JA1JejwiDtWWGdASahulcJMLOixZsiQeeOCBpNizZ88mb/fMmTPJsefOnUuOvXjxYnJs6nBOznZz/m5yJhtzho4uXbqUHJtj7ty5ybFjY2NJcWfPntXk5CSTmBhI15rE5F4oAFAoGjgAFIoGDgCFooEDQKFo4ABQKBo4ABSKBg4AhaKBA0ChaOAAUCgaOAAUqm/v77pu3brk2Jyx7JyFgl966aXk2KNHjybHpo7S92qMPWfx4Zw1MXu1CHRqDjm5AoOAI3AAKFSVNTFX2f4P26/ZPmT7C3UmBrSF2kYpqpxCmZL05xFx0PYiSQds742I12rKDWgLtY0idH0EHhEnIuJg5/F5Sa/rGusGAqWhtlGKWs6B214jabOkl+vYHtAvqG30s8pXodheKOlfJP1pRPzPNf58p6SdkrRgwYKqbwc0Jqe2gTZUOgK3PVczBf5kRHzjWjGzV+4eHR2t8nZAY3Jru9nsgBlVrkKxpK9Lej0i/rq+lIB2UdsoRZUj8F+R9AeSfsv2K53/fremvIA2UdsoQpVV6f9TUisLyAK9RG2jFI2O0o+OjmrNmjVJsTkrwi9cuDA5dtOmTcmxOaPsx48fT44dHx9Pisv5DHLG2Hs1cp6T7+TkZHJs6q0ScrYJDAJG6QGgUDRwACgUDRwACkUDB4BC0cABoFA0cAAoFA0cAApFAweAQtHAAaBQNHAAKFSjo/SLFi3S1q1bk2IPHDiQvN2DBw8mx+aMvK9bty45NvUWAZJ06NChpLhe3X733LlzybE5q8fnrHY/b9685NiJiYmkuJxcgUHAETgAFKpyA7c9ZPt7tv+1joSAfkFto9/VcQT+Bc0s+goMGmobfa3qkmorJT0g6Wv1pAP0B2obJah6BP4VSV+UlH4zaqAMXxG1jT5XZU3MT0o6GRE/83IR2ztt77e9P+fqB6At3dR2Q6kBP6Xqmpifsv0DSU9rZv3Af7w6aPbK3bfeemuFtwMak13bTScISBUaeER8KSJWRsQaSQ9K+veI+ExtmQEtobZRCq4DB4BC1TKJGRHflvTtOrYF9BNqG/2s0VH6kZERrV69Oil2aGgoebs5K7LnjNLnSL1FgCR98MEHSXE5X/rmjN2fP38+OTZnPD119XhJsp0cm1oLrEqPmw2nUACgUDRwACgUDRwACkUDB4BC0cABoFA0cAAoFA0cAApFAweAQtHAAaBQNHAAKFSjo/RDQ0NasGBBUuz4+Hjydu+///7k2MOHDyfHHjlyJDk251a5mzdvTorbt29f8jY//PDD5Nickffh4d6UyOXLl5NjU1ewz/kMgEHAETgAFKrqmpiLbT9j+7Dt123/Ul2JAW2itlGCqv8+/htJ34yI37M9ImmshpyAfkBto+913cBt3yrp1yR9VpIiYkLSRD1pAe2htlGKKqdQ1ko6JekfbH/P9tdsf+QbytkLv6beBxtoWXZtN58iUK2BD0v6RUl/FxGbJf2vpF1XB81e+HXJkiUV3g5oTHZtN50gIFVr4MckHYuIlzvPn9FM0QOlo7ZRhCqr0r8v6V3bGzovbZX0Wi1ZAS2itlGKqleh/LGkJzvf0r8t6Q+rpwT0BWobfa9SA4+IVyQln/+bM2eOxsbSrsbKWdQ4dbpTkhYvXpwce/r06eTYN954o/Yc7r777uRtHjx4MDk2ZxHoOXPS/5GWE5s6XSlJy5cvT4rLWaz5RnJre5DlLGzdCzkLYN9smMQEgELRwAGgUDRwACgUDRwACkUDB4BC0cABoFA0cAAoFA0cAApFAweAQtHAAaBQjS5qLKWP5eaM0ueMZecslpwTu2fPnuTYM2fOJMUtW7YseZurVq1Kjj127Fhy7OTkZHLswoULk2PvvPPO5Nh77rknKS7ndgYoR84o/802ds8ROAAUigYOAIWquir9n9k+ZPtV20/ZTj+XAfQxahsl6LqB214h6U8kbYmITZKGJD1YV2JAW6htlKLqKZRhSfNtD0sak/Re9ZSAvkBto+9VWVLtuKS/kvRDSScknYuIF66Om71yd84CCUBbuqntpnMEpGqnUG6TtEPSWknLJS2w/Zmr42av3L106dLuMwUa0k1tN50jIFU7hbJN0jsRcSoiJiV9Q9Iv15MW0CpqG0Wo0sB/KOl+22OeuXp+q6TX60kLaBW1jSJUOQf+sqRnJB2U9P3Oth6rKS+gNdQ2SlF1VfpHJT2aGm9bIyMjSbEXL15MziNn7H54OH2Xc0bZc8b5L1y4kBR39OjR5G2uXLkyOfa+++5Ljr3jjjuSY1evXp0cu3HjxuTYNWvWJMU9//zzydu8kdzaHmRtj6fnjNLfbGP3TGICQKFo4ABQKBo4ABSKBg4AhaKBA0ChaOAAUCgaOAAUigYOAIWigQNAoWjgAFCoRlelt508yn7lypXk7eaM0s+dOzc5dv369cmxW7ak31H07NmzSXGpq9dLUs691teuXZsce++99ybHbtu2LTk2ZwX71Nsq5NQBypEz8p4zSj8IOAIHgELdsIHbftz2SduvznrtY7b32n6z8/O23qYJ1I/aRulSjsB3S9p+1Wu7JL0YEeslvdh5DpRmt6htFOyGDTwiviPp6pOxOyQ90Xn8hKRP15sW0HvUNkrX7Tnw8Yg40Xn8vqTxmvIB2kZtoxiVv8SMma99r/vV7+yVu0+dOlX17YDG5NR2g2kBP9FtA/+R7WWS1Pl58nqBs1fuvv3227t8O6AxXdV2Y9kBs3TbwJ+T9Ejn8SOS9tSTDtA6ahvFSLmM8ClJ35W0wfYx25+T9GVJv237TUnbOs+BolDbKN0NxyIj4qHr/NHWmnMBGkVto3SNjtJPT0/r8uXLSbE5Y9HT09PJsRMTE8mxy5cvT459+OGHk2NXrFiRFPfOO+8kb3P+/Pk9iR0bG+vJdnNGnlP/fm+2MWp81CCsNJ+DUXoAKBQNHAAKRQMHgELRwAGgUDRwACgUDRwACkUDB4BC0cABoFA0cAAoFA0cAArV6Ch9RCSP0s+bNy95uzlj9znj1qOjo8mxd911V3Ls4sWLk+IuXLiQvM1Lly4lx05NTSXH5nwG58+fT47NGbtPHY++2caoAY7AAaBQ3a5K/5e2D9v+L9vP2l7c0yyBHqC2UbpuV6XfK2lTRPyCpDckfanmvIAm7Ba1jYJ1tSp9RLwQET8+kfqSpJU9yA3oKWobpavjHPgfSfq36/3h7IVfT58+XcPbAY1Jru0GcwJ+olIDt/0XkqYkPXm9mNkLvy5durTK2wGNya3t5jID/l/XlxHa/qykT0raGiyFggFCbaMUXTVw29slfVHSr0dE+sXKQJ+jtlGSblel/1tJiyTttf2K7b/vcZ5A7ahtlK7bVem/3oNcgEZR2yhd46P0qSuM54yR92KFcylvjDxn9D91lD41LlfOyPncuXOTYycmJpJjR0ZGkmPnzEn7rj3nlgrAIGCUHgAKRQMHgELRwAGgUDRwACgUDRwACkUDB4BC0cABoFA0cAAoFA0cAApFAweAQrnJu2XaPiXpv696eamkQV3pYVD3rV/36+cj4vY23vgmq+1B3S+pf/ftmrXdaAO/Ftv7B/WG+IO6b4O6X3Ub1M9pUPdLKm/fOIUCAIWigQNAofqhgT/WdgI9NKj7Nqj7VbdB/ZwGdb+kwvat9XPgAIDu9MMROACgC602cNvbbR+x/ZbtXW3mUifbP7D9/c6aivvbzqcK24/bPmn71Vmvfcz2Xttvdn7e1maO/Yja7n+DUNutNXDbQ5K+KukTkjZKesj2xrby6YHfjIh7S7ok6Tp2S9p+1Wu7JL0YEeslvdh5jg5quxi7VXhtt3kE/nFJb0XE2xExIelpSTtazAfXEBHfkXTmqpd3SHqi8/gJSZ9uMqcCUNsFGITabrOBr5D07qznxzqvDYKQ9ILtA7Z3tp1MD4xHxInO4/cljbeZTB+itstVVG03uir9TeRXI+K47Z+TtNf24c5v+4ETEWGbS5luHtR2H2nzCPy4pFWznq/svFa8iDje+XlS0rOa+Sf1IPmR7WWS1Pl5suV8+g21Xa6iarvNBr5P0nrba22PSHpQ0nMt5lML2wtsL/rxY0m/I+nVn/1/Fec5SY90Hj8iaU+LufQjartcRdV2a6dQImLK9uclfUvSkKTHI+JQW/nUaFzSs7almc/3nyLim+2m1D3bT0n6DUlLbR+T9KikL0v6Z9uf08wd+H6/vQz7D7VdhkGobSYxAaBQTGICQKFo4ABQKBo4ABSKBg4AhaKBA0ChaOAAUCgaOAAUigYOAIX6P0++YmDExNZZAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOX0lEQVR4nO3df4hd9ZnH8c+T+ZGZJJJMJnbya9IEjAUtxZRQ3a7sD80uyW5pCi6LQqv9Af2rPxFKSv/wX2FlaWHLLtK4EdYqxa1EFtoabEsprJIfDVvjWDXa1km1iTPGqJnfefrHXN3badJ8v/ece859bt4vkLn35vGc59x58szJmfs9j7m7AADxLKs7AQBAa2jgABAUDRwAgqKBA0BQNHAACIoGDgBB9Va5s76+Ph8YGEiKXbYs/WdLb2/6YSwsLCTHzs/Pt2W7qR/dvHDhQunbzI01s+TYduWbs013T0+4RGbG53HRVher7Uob+MDAgHbu3JkUOzg4mLzddevWJceePXs2OXZycrIt252enk6Km52dTd7m1NRUcmy7GvjMzExybOp7kCPnBy7QDQpdQjGz3Wb2KzN70cz2lZUUUDdqGxG03MDNrEfStyXtkXSdpDvM7LqyEgPqQm0jiiJn4B+R9KK7v+Tus5IekbS3nLSAWlHbCKFIA98k6ZWm5+ON1/6ImX3ezI6Y2ZG5ubkCuwMqk13blWUGNGn7xwjd/X533+nuO/v6+tq9O6AyzbVddy64MhVp4KckjTY939x4DYiO2kYIRRr4YUnbzWybmfVLul3S4+WkBdSK2kYILX8O3N3nzewLkn4kqUfSA+5+orTMgJpQ24jCqhzoMDo66nfffXdS7PHjx5O3m7MCcHh4ODm2v78/Obanpyc59qmnnkqKO3nyZPI2JyYmkmPbsYhGyls9m7NyNXOFJysx0ZUuVtvcCwUAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUJXOxFy9erX27NmTFDsyMpK83bGxseTYnNmRq1atSo69/vrrk2NTl7KfOpV+A7yc9yvnPchZxp6zPD5Har45xwV0A87AASCoIjMxR83sJ2b2rJmdMLMvl5kYUBdqG1EUuYQyL+ludz9mZldJOmpmh9z92ZJyA+pCbSOEls/A3f1Vdz/WePyWpDFdZG4gEA21jShKuQZuZlsl7ZD0dBnbAzoFtY1OVriBm9kqSf8t6Svufu4if/7e5O433nij6O6AyuTUdvXZAQUbuJn1abHAH3L3718spnly99DQUJHdAZXJre1qswMWFfkUiknaL2nM3f+1vJSAelHbiKLIGfhfSvqUpFvM7Hjjv38oKS+gTtQ2Qigylf7nkmoZIAu0E7WNKCpdSt/X16f169cnxd58883J2924cWNy7NGjR5NjjxxJ/93U+Ph4cuw111yTFLd169bkbZ44cSI5dvny5cmxOd58883kWPf0Ie450+6BKwl/MwAgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAEVelS+oWFBZ079ye3VS5sy5YtybE9PT3JsTkT2XMmyKe65ZZbkmMnJiaSY3OWvOcsu3/rrbeSY3OW0vf19SXFMZUeVxrOwAEgqDIm8vSY2S/M7H/KSAjoFNQ2Ol0ZZ+Bf1uLQV6DbUNvoaEVHqm2W9I+SvlNOOkBnoLYRQdEz8G9K+pqk9N/2ATF8U9Q2OlyRmZgfk3Ta3f/shITmyd2Tk5Ot7g6oTCu1XVFqwB8pOhPz42b2a0mPaHF+4H8tDWqe3L127doCuwMqk13bVScISAUauLt/3d03u/tWSbdL+rG7f7K0zICaUNuIgs+BA0BQpazEdPefSvppGdsCOgm1jU5W6VJ6STKzpLicJe/9/f3JsSMjI8mxN910U3Ls2Fj6x4Wff/75pLjVq1cnb3PHjh3JsYcPH06Offvtt5NjU5e8S1Jvb/mll1pbQLfgEgoABEUDB4CgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoCpdSr9s2TKtWLEiKTZnqfWyZek/h3K2u3LlyuTYNWvWJMemTpBPXXKfu/9rr702OfbYsWPJsRcupM8+yPmepcaylD4Od687ha6oF87AASCoojMx15jZo2b2nJmNmdlflJUYUCdqGxEUvYTyLUk/dPd/MrN+SWnXR4DOR22j47XcwM1staS/kvRpSXL3WUmz5aQF1IfaRhRFLqFsk3RG0n+a2S/M7Dtm9ie/9Wse/Jr6yzugZtm1XX2KQLEG3ivpw5L+3d13SHpH0r6lQc2DX4eHhwvsDqhMdm1XnSAgFWvg45LG3f3pxvNHtVj0QHTUNkIoMpX+NUmvmNkHGi/dKunZUrICakRtI4qin0L5oqSHGr+lf0nSZ4qnBHQEahsdr1ADd/fjkpKv//X29mpoaCgpdn5+PjmP6enp5NicYckDAwPJsTnDklNjDx48mLzNycnJ5NgNGzYkx46OjibHjo+PJ8fOzc0lx65atSopLmd15+Xk1jbiyVkN2qmrNlmJCQBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoSocau7tmZmaSYnOWsecM052amkqOzVmanbNEP3Upe857cP78+eTYkydPJsdu3rw5OfbGG29Mjl2/fn1y7JYtW5Li7rvvvuRtol6dsDQ9Zyl9py675wwcAIKigQNAUEWn0n/VzE6Y2TNm9rCZpf+bH+hg1DYiaLmBm9kmSV+StNPdPyipR9LtZSUG1IXaRhRFL6H0Sho0s15JKyT9rnhKQEegttHxioxUOyXpPkm/lfSqpDfd/Ymlcc2Tu8+cOdN6pkBFWqntqnMEpGKXUIYk7ZW0TdJGSSvN7JNL45ond1999dWtZwpUpJXarjpHQCp2CWWXpJfd/Yy7z0n6vqSPlpMWUCtqGyEUaeC/lXSTma2wxU+u3ypprJy0gFpR2wihyDXwpyU9KumYpF82tnV/SXkBtaG2EUXRqfT3SLonNd7M1N/fn7rt5DxylrEvLCwkx+YspU89Lknavn17UtzOnemXVs+ePZscmzPB/vXXX0+O3bZtW3LsDTfckBy7a9eupLj9+/cnb/Nycmsb8eQsec/pR1ViJSYABEUDB4CgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoCqdSn/hwgVNT08nxeYsTZ+dnU2OzVl2nzPtfmZmJjl248aNSXF33nln8jY3bdqUHPvyyy8nxw4ODrYldsWKFaVvN+fWB0COKifN56DiASCoyzZwM3vAzE6b2TNNr601s0Nm9kLj61B70wTKR20jupQz8AOSdi95bZ+kJ919u6QnG8+BaA6I2kZgl23g7v4zSUvvP7pX0oONxw9K+kS5aQHtR20julavgY+4+6uNx69JGikpH6Bu1DbCKPxLTF+80/kl73bePLk7ZzgAULec2q4wLeA9rTbw35vZBklqfD19qcDmyd3r1q1rcXdAZVqq7cqyA5q02sAfl3RX4/Fdkg6Wkw5QO2obYaR8jPBhSf8r6QNmNm5mn5N0r6S/M7MXJO1qPAdCobYR3WVXYrr7HZf4o1tLzgWoFLWN6CpdSm9m6uvrS4qdmppK3u78/Hxy7MDAQHJsu5bdp+aQM+X9tttuS449f/58cmzqrQ+kvO/D8uXLk2PPnTuXFLewsJC8TaAbsJQeAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCooEDQFCVLqVfWFhIXhadM2E8Zwn1O++8kxybIyeH1An2OcvN16xZ05bYHDmTu1NvqSBJs7OzSXG9vZWWM1A7zsABIKhWp9L/i5k9Z2b/Z2aPmdmatmYJtAG1jehanUp/SNIH3f1Dkp6X9PWS8wKqcEDUNgJraSq9uz/h7u/eO/QpSZvbkBvQVtQ2oivjGvhnJf3gUn/YPPh1YmKihN0BlUmu7QpzAt5TqIGb2TckzUt66FIxzYNfh4eHi+wOqExubVeXGfD/Wv7clZl9WtLHJN3q7l5aRkDNqG1E0VIDN7Pdkr4m6a/dPX0+F9DhqG1E0upU+n+TdJWkQ2Z23Mz+o815AqWjthFdq1Pp97chF6BS1Daiq3ztcer09pwp73Nzc8mxg4ODybE5ObRjKX3O9Pic9yDnNgU5BgYG2rLd1GPL+X4B3YCl9AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCsirvlmlmZyT9ZsnL6yS9XlkS1erWY+vU43q/u19dx46vsNru1uOSOvfYLlrblTbwizGzI916Q/xuPbZuPa6ydev71K3HJcU7Ni6hAEBQNHAACKoTGvj9dSfQRt16bN16XGXr1vepW49LCnZstV8DBwC0phPOwAEALai1gZvZbjP7lZm9aGb76sylTGb2azP7ZWOm4pG68ynCzB4ws9Nm9kzTa2vN7JCZvdD4OlRnjp2I2u583VDbtTVwM+uR9G1JeyRdJ+kOM7uurnza4G/d/YZIH0m6hAOSdi95bZ+kJ919u6QnG8/RQG2HcUDBa7vOM/CPSHrR3V9y91lJj0jaW2M+uAh3/5mkySUv75X0YOPxg5I+UWVOAVDbAXRDbdfZwDdJeqXp+XjjtW7gkp4ws6Nm9vm6k2mDEXd/tfH4NUkjdSbTgajtuELVduVT6a8QN7v7KTN7n6RDZvZc46d913F3NzM+ynTloLY7SJ1n4KckjTY939x4LTx3P9X4elrSY1r8J3U3+b2ZbZCkxtfTNefTaajtuELVdp0N/LCk7Wa2zcz6Jd0u6fEa8ymFma00s6vefSzp7yU98+f/r3Ael3RX4/Fdkg7WmEsnorbjClXbtV1Ccfd5M/uCpB9J6pH0gLufqCufEo1IeszMpMX397vu/sN6U2qdmT0s6W8krTOzcUn3SLpX0vfM7HNavAPfP9eXYeehtmPohtpmJSYABMVKTAAIigYOAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQfwBppnGdph3sPAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOP0lEQVR4nO3dfYxc5XXH8d/Z9Y537TUYe+ni19oIE8lEFa6sQArqC7jFSaw4ElUFUgppI+WvtGlVKXLUP/g3UqMolRq1Qgk1UimookFGlZpg0URRpID8EtRglgCGNKwLsfHKGNh3+/SPndDJxo6fZ+6de+8Zvh8J7cz4+Llndo8O13fvM8fcXQCAeAbqTgAA0B0aOAAERQMHgKBo4AAQFA0cAIKigQNAUCuqPFir1fKRkZGk2BUr0lO7cOFCcuzi4mJP1s25HfPixYulr5kTa2bJsam55uaQIzXfixcvyt3T31yJzIz7cdFTl6rtShv4yMiIbrvttqTYsbGx5HXPnTuXHDs1NdWTdWdnZ5Nj5+fnk+JmZmaS1+xVA5+bm0uOzfke5BgaGqr1+EBTFbqEYmZ7zewnZvaKmR0oKymgbtQ2Iui6gZvZoKSvS/qYpJ2S7jWznWUlBtSF2kYURc7APyLpFXd/1d3nJT0maX85aQG1orYRQpEGvknS6x3PJ9uv/RIz+5yZHTWzo6nXfoGaZdd2ZZkBHXp+G6G7P+juu919d6vV6vXhgMp01nbdueCDqUgDPyVpS8fzze3XgOiobYRQpIEfkbTDzLabWUvSPZKeLCctoFbUNkLo+j5wd180s89L+o6kQUkPufuJ0jIDakJtIwqrcqDD2NiY79u3Lyk2Zwfg+vXrk2NzrsMPDg4mxz7zzDPJsSdPnkyKO3v2bPKavdrEMjCQ/o+0nJ2rOT/f1F25i4uL7MRE37pUbfNZKAAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASAoGjgABEUDB4CgaOAAEBQNHACCqnQm5lVXXaW77rorKXZiYiJ53ZzZkaOjo8mxN910U3Jszlb2U6fSPthufHw8ec2c70HONvac7fE5cvJdWFhIisuZ9Qn0A87AASCoIjMxt5jZd83sBTM7YWZfKDMxoC7UNqIocgllUdLfuPtxM1sj6ZiZHXb3F0rKDagLtY0Quj4Dd/c33P14+/E7kiZ0ibmBQDTUNqIo5Rq4mW2TtEvSs2WsBzQFtY0mK9zAzWxU0r9L+it3P3+JP39/cvc777xT9HBAZXJqu/rsgIIN3MyGtFTgj7j7ty4V0zm5e82aNUUOB1Qmt7arzQ5YUuQuFJP0TUkT7v7V8lIC6kVtI4oiZ+C3SfpTSXeY2XPt/z5eUl5AnahthFBkKv0PJLH1DX2H2kYUlW6lHx0d1e23354Uu3HjxuR1jx07lhx79Gj675smJyeTY2+44Ybk2G3btiXFnThxInnNlStXJsfmePvtt5Nj3dMHs+dMux8eHk6KW1xcTF4T6AdspQeAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQVKVb6XNs3bo1OXZwcDA5Nmcie+r0+Fx33HFHUtzZs2eT18zZ8p6z7T7nM9xzttIPDQ0lx6ZOm5+enk5eE+gHnIEDQFBlTOQZNLMfmdl/lJEQ0BTUNpqujDPwL2hp6CvQb6htNFrRkWqbJX1C0jfKSQdoBmobERQ9A/+apC9KSv/NIBDD10Rto+GKzMTcJ+m0u//aaQqdk7unpqa6PRxQmW5qu6LUgF9SdCbmJ83sp5Ie09L8wH9ZHtQ5uXvdunUFDgdUJru2q04QkAo0cHf/krtvdvdtku6R9F/u/unSMgNqQm0jCu4DB4CgStmJ6e7fk/S9MtYCmoTaRpNVupXezLRiRdohW61W8rrj4+PJsbfeemty7MRE+i3AL730UnLs1VdfnRS3a9eu5DWPHDmSHPvuu+8mx+ZseU/92eaam5vrybpAdFxCAYCgaOAAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoGjgABBUpVvpBwYGtHr16uTYVDlbuFOPL0lr165Njs2ZIJ+67T7n+DfeeGNy7PHjx5NjL15Mn2eQ8zPLiR0eHk6KO3/+fPKa6A13r/X4Zlbr8avGGTgABFV0JuZaM3vczF40swkz+2hZiQF1orYRQdFLKH8v6dvu/sdm1pK0qoScgCagttF4XTdwM7ta0u9K+owkufu8pPly0gLqQ20jiiKXULZLOiPpn83sR2b2DTP7ld8Qdg5+zflFH1Cj7NquPkWgWANfIem3Jf2ju++S9J6kA8uDOge/rl+/vsDhgMpk13bVCQJSsQY+KWnS3Z9tP39cS0UPREdtI4QiU+nflPS6mX2o/dKdkl4oJSugRtQ2oih6F8pfSHqk/Vv6VyX9WfGUgEagttF4hRq4uz8nKfn638DAgEZGRpJiZ2dnk/MYHBxMjk3d1SflDUvOiT106FBS3NTUVPKaGzZsSI7dsmVLcuzk5GRy7MLCQnLs6Ohocuz111+fFHfmzJnkNa8kt7bRDDk7Qfth1yY7MQEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAEVelQYzPT0NBQUmzOMN2ZmZnk2Jxhujlb9HO2sqdu55+enk5e8+TJk8mxmzdvTo695ZZbkmOvu+665NitW7cmx+7cuTMpbmJiInlN9Ebd29NzttL3w7Z7zsABICgaOAAEVXQq/V+b2Qkze97MHjWz9I/6AxqM2kYEXTdwM9sk6S8l7Xb3D0salHRPWYkBdaG2EUXRSygrJI2Y2QpJqyT9b/GUgEagttF4RUaqnZL0FUk/k/SGpLfd/anlcZ2Tu8v8wH2gV7qp7apzBKRil1CukbRf0nZJGyWtNrNPL4/rnNx97bXXdp8pUJFuarvqHAGp2CWUPZJec/cz7r4g6VuSfqectIBaUdsIoUgD/5mkW81slS3d5X6nJHZSoB9Q2wihyDXwZyU9Lum4pB+313qwpLyA2lDbiKLoVPoHJD2Q+XeS4nK2sV+4cCE5NmcrfavVSo7dsWNHcuzu3WmXTM+dO5e8Zs4E+7feeis5dvv27cmxN998c3Lsnj17kmNTJ9jnTLq/km5qG/XL2fKes5W+qdiJCQBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoSqfSu3vytvf5+fnkdXO23edMu5+bm0uO3bhxY3LsfffdlxS3adOm5DVfe+215NiRkZGexK5ataon6/bDlmc0T1MnzefgDBwAgrpiAzezh8zstJk93/HaOjM7bGYvt79e09s0gfJR24gu5Qz8oKS9y147IOlpd98h6en2cyCag6K2EdgVG7i7f1/S8s8q3S/p4fbjhyV9qty0gN6jthFdt9fAx939jfbjNyWNl5QPUDdqG2EU/iWmL90icNnbBDond+cMEgDqllPbFaYFvK/bBv5zM9sgSe2vpy8X2Dm5e2xsrMvDAZXpqrYryw7o0G0Df1LS/e3H90s6VE46QO2obYSRchvho5J+KOlDZjZpZp+V9GVJf2hmL0va034OhEJtI7or7sR093sv80d3lpwLUClqG9FVvpV+ZmYmKXZxcTF53eHh4eTYXm27z8khddL73Xffnbzm9PR0cuzs7GxybM7PYeXKlcmx58+fT45N3aKf8/MC+gFb6QEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAEVflW+tSt2anT6yXpvffe6zal0nLImWCfuuV87dq1yWvmxObImdw9NDSUHDs/P1/6ujkfkwD0A87AASCobqfS/52ZvWhm/21mT5jZ2p5mCfQAtY3oup1Kf1jSh939tyS9JOlLJecFVOGgqG0E1tVUend/yt1/cTH7GUmbe5Ab0FPUNqIr4xr4n0v6z8v9IUONEVhybVeYE/C+Qg3czP5W0qKkRy4Xw1BjRJRb29VlBvy/rm8jNLPPSNon6U5399IyAmpGbSOKrhq4me2V9EVJv+fu6bO8gIajthFJt1Pp/0HSGkmHzew5M/unHucJlI7aRnTdTqX/Zg9yASpFbSO6xm6lX1hYSF53ZGQkOTZncnmvttKnTpDP+R4MDPRmU+3w8HBP1s15b61Wqyc5ANGxlR4AgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUFblp2Wa2RlJ/7Ps5TFJ/TrpoV/fW1Pf12+6+7V1HPgDVtv9+r6k5r63S9Z2pQ38UszsaL9+IH6/vrd+fV9l69fvU7++Lynee+MSCgAERQMHgKCa0MAfrDuBHurX99av76ts/fp96tf3JQV7b7VfAwcAdKcJZ+AAgC7U2sDNbK+Z/cTMXjGzA3XmUiYz+6mZ/bg9U/Fo3fkUYWYPmdlpM3u+47V1ZnbYzF5uf72mzhybiNpuvn6o7doauJkNSvq6pI9J2inpXjPbWVc+PfAH7n5zpFuSLuOgpL3LXjsg6Wl33yHp6fZztFHbYRxU8Nqu8wz8I5JecfdX3X1e0mOS9teYDy7B3b8vaWrZy/slPdx+/LCkT1WZUwDUdgD9UNt1NvBNkl7veD7Zfq0fuKSnzOyYmX2u7mR6YNzd32g/flPSeJ3JNBC1HVeo2q50Kv0HyO3ufsrMfkPSYTN7sf1/+77j7m5m3Mr0wUFtN0idZ+CnJG3peL65/Vp47n6q/fW0pCe09E/qfvJzM9sgSe2vp2vOp2mo7bhC1XadDfyIpB1mtt3MWpLukfRkjfmUwsxWm9maXzyW9EeSnv/1fyucJyXd3358v6RDNebSRNR2XKFqu7ZLKO6+aGafl/QdSYOSHnL3E3XlU6JxSU+YmbT0/f1Xd/92vSl1z8welfT7ksbMbFLSA5K+LOnfzOyzWvoEvj+pL8PmobZj6IfaZicmAATFTkwACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUP8HUYhxVa5/I9QAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAC3CAYAAAALgwWHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOcElEQVR4nO3dbYxc9XXH8d/ZZ3Zt2dhLFz/WBkwkY0W4sgJ9UJ/sVk6J6ryoKpBSSBrJr9KmVaXIUV/wNlKrKpUatUIJNVIpqKJBhkpJsGijqFJBtjeowWADhiSsC7G9KztLbe/sw+mLnaQTY8f//9w7986Z/X4ktDPjw73nLkfHlzv33L+5uwAA8fTVnQAAoD00cAAIigYOAEHRwAEgKBo4AARFAweAoAaq3NnIyIiPjY0lxS4uLiZvd2FhITk2Z7s5t1guLS2Vvt2c/ZtZcmwncs2Vk29fX9p5xuLiopaWltI3XCIz435cdJS7f6i2K23gY2NjeuCBB5JiL168mLzdmZmZ5NhLly4lx165ciU5ttFolL7dTjXwubm55NirV68mx+YYHBxMjh0dHU2Ky6kZoBcUuoRiZvvN7LSZvWVmh8pKCqgbtY0I2m7gZtYv6SuSPi5pp6SHzGxnWYkBdaG2EUWRM/CPSXrL3d9294akpyUdKCctoFbUNkIo0sA3SXq35f1U87OfYWYHzey4mR3PufYK1Ci7tivLDGjR8dsI3f0xd9/j7nuGh4c7vTugMq21XXcuWJmKNPCzkra0vN/c/AyIjtpGCEUa+DFJO8xsu5kNSXpQ0nPlpAXUitpGCG3fB+7uC2b2OUnfktQv6XF3P1laZkBNqG1EYVUu6LB+/XpPHeQZHx9P3m7OUEh/f39y7EsvvZQce+bMmeTY6enppLhODdGkTjZKeZOrOROeAwPp5w6p/80ajQaTmOhZ15vE5FkoABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKqdE3M4eFhbdu2LSk2Zz3KVatWJcfu2rUrOTZnlP3s2fSH1U1MTCTF5fwOcsbYc8bjc+TkOz8/nxyb+qiEnG0CvYAzcAAIqsiamFvM7D/M7DUzO2lmny8zMaAu1DaiKHIJZUHSX7j7pJmtlnTCzI66+2sl5QbUhdpGCG2fgbv7e+4+2Xw9K+l1XWfdQCAaahtRlHIN3My2Sdot6eUytgd0C2ob3azwXShmtkrSv0r6M3f/8XX+/KCkg5K0Zs2aorsDKpNT20AdCp2Bm9mglgv8SXf/+vViWlfuHh0dLbI7oDK5tV1tdsCyInehmKSvSXrd3f+mvJSAelHbiKLIGfivSvojSb9tZq80//m9kvIC6kRtI4Qiq9L/p6RaFpAFOonaRhSVjtKvXr1ae/fuTYo9ceJE8nYnJyeTY3NG3u+6667k2NRHBEjSyZMnk+KGh4eTt5nj0qVLybHu6Yut56x2PzIykhzbaDSS4nJyBXoBo/QAEBQNHACCooEDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqp0lH5oaEhbt25Niu3v70/ebs6K7Dmj9DlSHxEgSdPT00lxOSPvOWP3s7OzybE54+mpq8dL0vID/9Kk1gKr0mOl4QwcAIIq3MDNrN/Mvmtm/1ZGQkC3oLbR7co4A/+8lhd9BXoNtY2uVnRJtc2SHpD01XLSAboDtY0Iip6Bf1nSFySlf4sIxPBlUdvockXWxPyEpHPu/nNXXjCzg2Z23MyOz8zMtLs7oDLt1HZFqQE/o+iamL9vZt+X9LSW1w/8p2uDWlfuXrduXYHdAZXJru2qEwSkAg3c3b/o7pvdfZukByX9u7t/qrTMgJpQ24iC+8ABIKhSJjHd/duSvl3GtoBuQm2jm1U6St/f36+xsbGk2ImJieTt3n///cmxp06dSo49ffp0cuyaNWuSY3fv3p0Ud+zYseRtfvDBB8mxOSPvAwOdKZG5ubnk2NQV7HN+B0Av4BIKAARFAweAoGjgABAUDRwAgqKBA0BQNHAACIoGDgBB0cABICgaOAAERQMHgKAqHaXv6+vT6OhoUmzOqvSp4/mStHbt2uTYCxcuJMe+8cYbpedw9913J29zcnIyOXZpKX2Ngr6+9L/jc2JTx+MlaePGjUlxs7OzydtEOnevdf9mVuv+uxln4AAQVNE1Mdea2TNmdsrMXjezXy4rMaBO1DYiKHoJ5W8lfdPd/8DMhiSlXR8Buh+1ja7XdgM3szWSfl3SpyXJ3RuSGuWkBdSH2kYURS6hbJd0XtI/mtl3zeyrZvahbxNbF37N+VIQqFF2bVefIlCsgQ9I+iVJf+/uuyX9r6RD1wa1Lvw6Pj5eYHdAZbJru+oEAalYA5+SNOXuLzffP6Plogeio7YRQpFV6d+X9K6ZfaT50V5Jr5WSFVAjahtRFL0L5U8kPdn8lv5tSZ8pnhLQFahtdL1CDdzdX5GUdf0vdaorZxIzZ6ovZ7HknNgjR44kx87MzCTFbdiwIXmbW7ZsSY6dmppKjp2fn0+OXbVqVXLsHXfckRx7zz33JMXlTMPeTDu1jc7ImQRdaVObTGICQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACKrSRY3NTENDQ0mxV65cSd5uztj9wED6IeeMsueM81++fDkp7syZM8nb3Lx5c3Lsfffdlxx7++23J8du3bo1OXbnzp3Jsdu2bUuKe/7555O3iXR1j6fnjNKvtLF7zsABICgaOAAEVXRV+j83s5Nm9qqZPWVm6dcRgC5GbSOCthu4mW2S9KeS9rj7Lkn9kh4sKzGgLtQ2oih6CWVA0i1mNiBpVNL/FE8J6ArUNrpekSXVzkr6a0k/lPSepEvu/sK1ca0rd58/f779TIGKtFPbVecISMUuodwq6YCk7ZI2Shozs09dG9e6cvdtt93WfqZARdqp7apzBKRil1D2SXrH3c+7+7ykr0v6lXLSAmpFbSOEIg38h5LuN7NRW74jfq+k18tJC6gVtY0QilwDf1nSM5ImJX2vua3HSsoLqA21jSiKrkr/qKRHU+PNLHmUfXFxMTmPnFH6wcHB5NgdO3Ykx+7Zk34Z9OLFi0lxqavXS9KFCxeSY7dv354ce++99ybH7tu3Lzk2ZwX71Mcq5NTBzeTWNjonZ+Q9Z5S+FzCJCQBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoSlelX1pa0tzcXFJszlj00tJScmyj0UiO3bhxY3Lsww8/nBy7adOmpLh33nkneZu33HJLR2JHR0c7st2ckefU/74rbYwaH9YLK83n4AwcAIK6aQM3s8fN7JyZvdry2TozO2pmbzZ/3trZNIHyUduILuUM/LCk/dd8dkjSi+6+Q9KLzfdANIdFbSOwmzZwd/+OpGufa3pA0hPN109I+mS5aQGdR20junavgU+4+3vN1+9LmigpH6Bu1DbCKPwlpi9/9X/Dr/9bV+7OWXQAqFtObVeYFvBT7TbwH5nZBklq/jx3o8DWlbvHx8fb3B1QmbZqu7LsgBbtNvDnJD3SfP2IpCPlpAPUjtpGGCm3ET4l6b8kfcTMpszss5K+JOl3zOxNSfua74FQqG1Ed9NJTHd/6AZ/tLfkXIBKUduIrtJRendPHqUfGRlJ3m7O2H3OuPXw8HBy7J133pkcu3bt2qS4y5cvJ2/z6tWrybELCwvJsTm/g9nZ2eTYnLH71PHolTZGDTBKDwBB0cABICgaOAAERQMHgKBo4AAQFA0cAIKigQNAUDRwAAiKBg4AQdHAASCoykfpU1cYzxkj78QK51LeGHnO6H/qKH1qXK6ckfPBwcHk2EajkRw7NDSUHNvXl3aekfNIBaAXcAYOAEG1uyr9X5nZKTP7bzN71szWdjRLoAOobUTX7qr0RyXtcvePSnpD0hdLzguowmFR2wisrVXp3f0Fd//JM0lfkrS5A7kBHUVtI7oyroH/saRv3OgPWxd+nZ6eLmF3QGWSa7vCnICfKtTAzewvJS1IevJGMa0Lv65fv77I7oDK5NZ2dZkB/6/t2wjN7NOSPiFpr+fcxwd0OWobUbTVwM1sv6QvSPoNd0+/YRvoctQ2Iml3Vfq/k7Ra0lEze8XM/qHDeQKlo7YRXbur0n+tA7kAlaK2EV3Xrko/OjqavN3FxcXk2JwV2XPkjP7Pz88nxaWOkOfKWRE+R+pxSXmPKUgd52dVeqw0jNIDQFA0cAAIigYOAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACMqqfFqmmZ2X9INrPh6XdKGyJKrVq8fWrcf1i+5+Wx07XmG13avHJXXvsV23titt4NdjZsd79YH4vXpsvXpcZevV31OvHpcU79i4hAIAQdHAASCobmjgj9WdQAf16rH16nGVrVd/T716XFKwY6v9GjgAoD3dcAYOAGhDrQ3czPab2Wkze8vMDtWZS5nM7Ptm9r3mmorH686nCDN73MzOmdmrLZ+tM7OjZvZm8+etdebYjajt7tcLtV1bAzezfklfkfRxSTslPWRmO+vKpwN+y93vjXRL0g0clrT/ms8OSXrR3XdIerH5Hk3UdhiHFby26zwD/5ikt9z9bXdvSHpa0oEa88F1uPt3JM1c8/EBSU80Xz8h6ZNV5hQAtR1AL9R2nQ18k6R3W95PNT/rBS7pBTM7YWYH606mAybc/b3m6/clTdSZTBeituMKVduVrkq/gvyau581s1+QdNTMTjX/tu857u5mxq1MKwe13UXqPAM/K2lLy/vNzc/Cc/ezzZ/nJD2r5f+l7iU/MrMNktT8ea7mfLoNtR1XqNqus4Efk7TDzLab2ZCkByU9V2M+pTCzMTNb/ZPXkn5X0qs//98K5zlJjzRfPyLpSI25dCNqO65QtV3bJRR3XzCzz0n6lqR+SY+7+8m68inRhKRnzUxa/v3+s7t/s96U2mdmT0n6TUnjZjYl6VFJX5L0L2b2WS0/ge8P68uw+1DbMfRCbTOJCQBBMYkJAEHRwAEgKBo4AARFAweAoGjgABAUDRwAgqKBA0BQNHAACOr/AHM3YEHh0M8HAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"n = 0\n",
"for x,y,z in zip(*np.where(fiss==1)):\n",
" print(x,y,z)\n",
" s = (x-7,y-7,z-7)\n",
" e = (x+7,y+7,z+7)\n",
" try:\n",
" cube_img = img[s[0]:e[0],s[1]:e[1],s[2]:e[2]]\n",
" cube_fiss = fiss[s[0]:e[0],s[1]:e[1],s[2]:e[2]]\n",
" plt.figure(n)\n",
" plt.subplot(121)\n",
" plt.imshow(cube_img[7,:,:],cmap='gray')\n",
" plt.subplot(122)\n",
" plt.imshow(cube_fiss[7,:,:],cmap='gray')\n",
" except:\n",
" pass\n",
" n+=1\n",
" if n > 6:\n",
" break"
]
},
{
"cell_type": "code",
"execution_count": 135,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
echo 'pip intalling stuff'
pip install -r requirements.txt
echo 'downloading 2d image'
curl --output ok.png https://i.imgur.com/euaZxKf.png
echo 'downloading 3d image'
python utils.py
pydicom
SimpleITK
imageio
scipy
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
import os
import sys
import requests
import json
import zipfile
import pydicom
import SimpleITK as sitk
def download_images(series_instance_uid,image_root="."):
url = f"https://services.cancerimagingarchive.net/services/v4/TCIA/query/getImage?SeriesInstanceUID={series_instance_uid}"
zip_file_path = os.path.join(image_root,series_instance_uid+'.zip')
r = requests.get(url, allow_redirects=True)
print(r.status_code)
if r.status_code != 200:
raise LookupError(f"ohoh {r.status_code}!")
open(zip_file_path, 'wb').write(r.content)
folder_path = os.path.join(image_root,series_instance_uid)
with zipfile.ZipFile(zip_file_path, 'r') as zip_ref:
zip_ref.extractall(folder_path)
return folder_path
def imread(folder_path):
file_list = [os.path.join(folder_path,x) for x in os.listdir(folder_path)]
dicom_list = []
for image_file in file_list:
ds=pydicom.dcmread(image_file)
dicom_list.append((ds.InstanceNumber,image_file))
dicom_list = sorted(dicom_list,key=lambda x:x[0])
dicom_names = [x[1] for x in dicom_list]
reader = sitk.ImageSeriesReader()
reader.SetFileNames(dicom_names)
img_obj = reader.Execute()
return img_obj
if __name__ == "__main__":
series_instance_uid='1.3.6.1.4.1.14519.5.2.1.6279.6001.113679818447732724990336702075'
image_root = "/tmp/myimages"
os.makedirs(image_root,exist_ok=True)
folder_path = download_images(series_instance_uid,image_root=image_root)
img_obj = imread(folder_path)
writer = sitk.ImageFileWriter()
writer.SetFileName('ok.nii.gz')
writer.SetUseCompression(True)
writer.Execute(img_obj)
@pangyuteng
Copy link
Author

pangyuteng commented Mar 3, 2022

    ksizes = (1,100,100,1)
    strides = (1,5,5,1)
    rates = [1,1,1,1]
    padding = "SAME"
    
    input_shape = (1,100,100,1)
    volume = np.random.random(input_shape)
    print(volume.shape)
    i = tf.image.extract_patches(
        volume, ksizes, strides, rates, padding, name=None
    )


    print('input_shape',input_shape)
    print('ksizes',ksizes)
    print('strides',strides)
    print('input_shape[1]*input_shape[2]',float(input_shape[1])*float(input_shape[2]))
    print('output_shape',i.shape)
    
input_shape (1, 100, 100, 1)
ksizes (1, 100, 100, 1)
strides (1, 5, 5, 1)
input_shape[1]*input_shape[2] 10000.0
output_shape (1, 20, 20, 10000)
    

@pangyuteng
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment