quantitative_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/xiangtiange/opt/miniconda3/envs/chris/lib/python3.6/site-packages/dicom/__init__.py:53: UserWarning: \n",
      "This code is using an older version of pydicom, which is no longer \n",
      "maintained as of Jan 2017.  You can access the new pydicom features and API \n",
      "by installing `pydicom` from PyPI.\n",
      "See 'Transitioning to pydicom 1.x' section at pydicom.readthedocs.org \n",
      "for more information.\n",
      "\n",
      "  warnings.warn(msg)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from dipy.core.gradients import gradient_table\n",
    "from dipy.data import get_fnames\n",
    "from dipy.io.gradients import read_bvals_bvecs\n",
    "from dipy.io.image import load_nifti, save_nifti\n",
    "from dipy.segment.mask import median_otsu\n",
    "from dipy.reconst.dti import TensorModel\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = '...' # **YOUR ROOT OF DENOISED DATA**\n",
    "denoised_data = {}\n",
    "\n",
    "hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\n",
    "data, affine = load_nifti(hardi_fname)\n",
    "raw_max =  np.max(data, axis=(0,1,2), keepdims=True)\n",
    "bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\n",
    "\n",
    "# non b0 mask\n",
    "sel_b = bvals != 0\n",
    "b0_size = 10\n",
    "\n",
    "# create a new gradient table on non b0 directions\n",
    "gtab = gradient_table(bvals[sel_b], bvecs[sel_b])\n",
    "raw_gtab = gradient_table(bvals, bvecs)\n",
    "\n",
    "# PLEASE DOUBLE CHECK THE PATH BELOW\n",
    "# THE DENOISED DATA MUST CONTAIN B0 DIRETCIONS!!\n",
    "denoised_data['n2s'], _ = load_nifti(root + 'n2s_denoised.nii.gz')\n",
    "denoised_data['n2n'], _ = load_nifti(root + 'n2n_denoised.nii.gz')\n",
    "denoised_data['dip'], _ = load_nifti(root + 'dip_denoised.nii.gz')\n",
    "denoised_data['p2s'], _ = load_nifti(root + 'p2s_denoised.nii.gz')\n",
    "denoised_data['our'], _ = load_nifti(root + 'ours_denoised.nii.gz')\n",
    "\n",
    "# normalize denoised data\n",
    "for k in denoised_data.keys():\n",
    "    denoised_data[k] -= np.min(denoised_data[k], axis=(0,1), keepdims=True)\n",
    "    denoised_data[k] /= np.max(denoised_data[k], axis=(0,1), keepdims=True)\n",
    "    denoised_data[k] = np.concatenate((data[:,:,:,:b0_size], denoised_data[k]), axis=-1)\n",
    "\n",
    "print('Data loading done!')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing brain mask...\n",
      "Computing tensors...\n"
     ]
    }
   ],
   "source": [
    "print('Computing brain mask...')\n",
    "b0_mask, mask = median_otsu(data, vol_idx=[0])\n",
    "\n",
    "print('Computing tensors...')\n",
    "tenmodel = TensorModel(raw_gtab)\n",
    "tensorfit = tenmodel.fit(data, mask=mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing worst-case/best-case SNR using the corpus callosum...\n",
      "(81, 106, 76)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fe0b17530b8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAABpCAYAAACZFPC4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABAL0lEQVR4nO29ebxdZX3v/37WvOd95jEhgYREZhEFnKAlDr1KFWmLgFpbRflpq5bb+rr6u1pvW731+vJXtQ6lVK9TBXGmzpUC1oIgIoohiQnJSXLmc/Y8r+n5/bH2Wh6GHDKcgeh6v1555WSvtZ/1POvsfPZ3fZ/vIKSUxMTExMSsDcp6TyAmJibmt4lYdGNiYmLWkFh0Y2JiYtaQWHRjYmJi1pBYdGNiYmLWkFh0Y2JiYtaQWHR/SxFCfFoI8XfrPY+YmN821k10hRDXCCHuF0LUhRAzQojvCCGe2z12uhDiS0KIRSFERQjxCyHEDUIIdb3muxzLraV7/KRaT0xMzOqxLqIrhLgB+BDwPmAI2Ah8HHiZEOI04F7gMHC2lDIH/CFwAZBZj/kux3Jr6R4/qdYTExOzykgp1/QPkAPqwB8e4fjngW+t9bxWYy3Hsx7gUmASeDswD8wALwf+G/AroAi8c8n5zwLuAcrdcz8KGN1jAviH7jhV4CHgrO6xTwN/1/05A9wBfAQQ631f4z/xn9/kP+th6V4MWMDXjnB8B/DltZvOCfFka4HjW89wd9wx4N3ATcCrgGcAzwPeJYTY3D3XA/4C6O/O5zLgTd1jLwSeD5xO8AXxR0Bh6YWEEH3A7cB/SSnfIqWM88JjYlaR9RDdPmBRSukuc3xmDedzIjzZWsJzjnU9DvBeKaUD3EIgqB+WUtaklDuBh4FzAaSUP5VS/lhK6UopJ4AbgUuWjJMBthNYsLuklEvnMgrcBXxJSvk/j3GOMTExx8F6iG4B6BdCaMscH1nD+ZwIT7aW8JxjXU9BSul1f251/55bcrwFpCHapPumEGJWCFEl8C33A0gp/4PA3fAxYF4I8c9CiOyScV4CJIB/Osb5xcTEHCfrIbr3AB0CP+UT8QPgyjWbzYnxZGuB1V/PJ4DdwFYpZRZ4J4EvFwAp5UeklM8AziBwM/zVkvfeBHwX+LYQIrWKc4yJiemy5qIrpawQ+Ck/JoR4uRAiKYTQhRC/J4T4P8BfA88WQnxACDEMIITYIoT4vBAiv9bzXY6jWAus/noyBJtkdSHEduD/CQ8IIZ4phLhQCKEDDaAN+I95/58Be4B/E0IkVmA+MTExy7AuIWNSyg8CNwD/E1ggCKf6M+DrUspHCDaENgE7hRAV4CvA/UBtPea7HMutpXt8tdfzl8A13bFuAr645Fi2+1oJOEjg6vjAY+YvgTcQREx8QwhhrcCcYmJijoCIN6tjYmJi1o44DTgmJiZmDYlFNyYmJmYNiUU3JiYmZg2JRTcmJiZmDVkuqJ8XKH/4W7HL9u/+l8STnxUghPituCdSyqO+JzExMUdPbOnGrCnd8penrvc8jhUhxGuFED9a73msBEKIS4UQk+s9j5XiZKsNHYtuzJoipUxLKfev9zxifvsQQtwphHj9es8jFt2YmJiYNSQW3ZhjRggxIYT4y24HjIoQ4otLM9mEENcJIfYJIYpCiNuEEKNLjkkhxJbuz/9NCPGwEKImhJgSQvzlkvNeKoR4UAhRFkLcLYQ4Z5n5SCHEm4QQe7tj/a0Q4rTu+6pCiFuFEEb33J5ugaAFIUSp+/P4krFeK4TY3x3ngBDi2iNc8wNCiB8JIXIndjefcOwJIcRfde9vQwjxSSHEULcjSU0I8QMhRE/33C91ix1VhBA/FEKcuWScI97fx1zvLd3zxp/o+FMNIcTThRAPdNf1RYIyqMv+boUQ7yUoi/rRrovro93XPyyEONz9nPxUCPG81Z5/LLoxx8sfAS8GNgPnAK8FEEL8LvC/u8dHCNKPbznCGJ8E3iilzABnAf/RHePpwKeANxKUxrwRuE0IYS4znxcR1Bu+iKAA/D8T1CDe0B376u55CvB/gVMIuny0CCqxIYKiPx8Bfq87p2cDDy69iBBCEULc1F3zC7v1N1aDK4EXEBQpuhz4DkExo4HuGt7SPe87wFZgEHgA+NclYzzh/V2KEOLdBL+7S6SUT3k/b/fL8+vA54Be4Ev8uqDUEX+3Usr/F/hP4M+6Lq4/677nJ8B53bG+AHxptVPhY9GNOV4+IqWcllIWgX8j+OACXAt8Skr5gJSyA7wDuFgIsekJxnCAM4QQWSllSUr5QPf1NwA3SinvlVJ6UsrPEFRzu2iZ+fwfKWW1W2/4l8D3pZT7u6L4HeDpAFLKgpTyK1LKppSyBryXX9cfhqAg0FlCiISUcqY7XogO3EzwH/RyKWXzqO7U8fGPUso5KeUUgVjcK6X8mZSyTVA0P1zPp7p1ljvAe4Bzl1jfR7q/AEII8f8RFLr/HSnlwiquZSW5iOD38CEppSOl/DKBcB7N7/ZxSCk/332f262jYgLbVnMB6yu6Qjz6Z0UN/l7ys2JZKMkkSjIZvBbzVGF2yc9NuvV9CQqjHwwPSCnrBIV2xp5gjCsJ2hAdFELcJYS4uPv6KcB/77oWykKIMoHFOvoEY4Q8tt7wkeoPJ4UQNwohDoqg/vAPgbwQQpVSNoCrgOuBGSHEt0RQuS1kC0Hvu/8lpbSXmctK8KTrEUKoQoi/F0I80l3LRPd4f/fvI91fgDzBl9v/XkVrfTUYBaYe0+HkICz/uz3SYCJwk+3qumfKBB1W+o90/kqw9qLbFVphmijJJELTIlHVRodR+3pRMxnUvl6EpiOlBCHw2x2EqgbvS6UeLdgxTyWmCUQTiB7Z+4Cpx54opfyJlPJlBI/GXwdu7R46TNA5I7/kT1JKefMKzO+/E1gyF3brDz8/nGp3Tt+TUr6AwDWym6BKW8gu4E+A7wghVtUaOkquIfgS2EEgFpu6r4drOdL9haDy3EuB/yuEeM4azXclmAHGhHiUAGzs/r3s7xZ4VIx913/7dgJXWI+UMg9Ulpy/KqyZ6ApNCwTTNAOL1fPwm83g9VQKFAVZrYHtIF0Xb3ER6XkIIUBKhKqC9BFCIF0XJZFAaMvmdsSsDzcDfyKEOK/rg30fwaPxxNKThBCGEOJaIUSu25aoyq9r/d4EXC+CWsBCCJESQrxECLES3ZMzBJZiWQjRS1DvOJzTkBDiZd0vig5B09FH1R/uCv87gR+IoNPzepIhmGcBSBLca+BJ7y8AUso7CdxBXxVCPGvNZn1i3AO4wFtEULv6FQTNWWGZ322XOWBpjHimO9YCoHX921lWmdUX3a6rQGgaeF7QDdMyQVVRe3sQ4yOITAolncKrN/DqjcCq1fTg7YkE+D5KKhBZ6fmBxWvoCMOIXQ5PMaSUPwDeRVAzeAY4DXjlEU5/NTDRfRS8nkAAkFLeD1xHsAlSAvbR3ahbAT5E0KJoEfgxQeeMEIWgNvI0QdflS1hSFD6k62P+G+A/juCrXis+S/BoPUXQN+/Hjzn+hPd3KVLKfwf+lKCI/fmrO90Tp+vWeQXB56FI4A76avfwhzjy7xbgw8AfdCMbPgJ8r3vOrwjuY5vgKWtVWbae7kqkASuWFQiloaPkc2Aa+AsFcBzYsgmlXMOvVBGWhV+uIF0HxTQDt4IvUU4Zg3IVv1YPBvQlwjIRhg5CIOsNUBT85vHvacRpwI8nTgOOiVkdVvf5XAjQddT+LNIy8C0TFFCcHLJWQx44jC8lfquFsG2EoSNdBwC1tweZz1Db3kt6r4lotpCeh3LqRrxcAuH4qDOLeLaDkrAQponsdFZ1OTExMTEnyuq4F7ruBGEYCMvCG8jT3NaP8DzE1DyyVkP05BGqipLNoGYygftBWTIdVcVPGniGoLY9hxgdQslmaWztpTmWRPg+0nEQisBvtcHzEOZyYZwxMTEx68+qWbrCMMD3ke02yuFZLHWE9oYciY6Dn7Twkjpqs4Xsy1N79iaSk02Ujouy5wBKTx6ZSyMcj+wjdeqbUnj9GVQh0Joe1kQBf6EQbLQZBth2dD1iizcmJuYpzMqLrhCouSyiJwcdG9npIJstlHKDRLGOXCzin7EJX1MQ4wMoTRtrvoPa6IY9ZjLYpw3jpjTMYgdcH18TKG0XimUsx8UvlAKfsKIEFrIQyE4n2qyLiYmJeaqysqLbjSTwm03UZAK/J4tSawT+Vl2DdiCsaqWF35eiPZgkOWGj750GX4Jjg6qi1W30YhNRbeD3ZRE+SE0JwskKJcRQPxg6UlEQzXbg79U1UNUgxEwoyFoNv91e0eXFLMsxbTBKKfn+97/Pe97zHn7848duuq8cuq6zYcMGXvSiF3Hdddfx9Kc//USHjDddH8OxbLoe6Z4IIRBCoOs6QggURcFxHDRNQ1VVXNdFVVUSiQSKoqCqKp1OB8cJ9oAURUEIgaZpGIaB4zjYto2u63ieR6fTQdeDuH/XdUmlUrRaLRzHid7neR6madJsNlFVFd/30TQNKSWqqtJqtfB9H9/3n2gJR31PVjR6QZgmSiaNbLWD2NpkAjk8gNJogZTIVhthGMhMEjefwDdUjPk6olTFr9YQVjeG13UDy9U0ob+H+pn9qC2fxEQZihXoz1N8Rh/Cg8yhNvpMObCqmy1EKom3sIgwDPza0XU4j6MXHs9xRC8c9X3xfZ+PfOQjvOMd76C9hl+MQ0ND/P7v/z5//ud/ztlnn328w8SflcdwoqKrdePthRCYponrupGISikxTRPHcdB1HV3XyWQyOI6DYRh0Oh1c18XrxvSn02mklLRaLVRVxbZthBA0Go3oWCaTodPp0Gg0SCQSeJ5Hu93GMAw8z8O27UhcpZRIKUkkEti2jeM4jzoWzvuxOrr6oisESiIRWJobR/ENDbVUo7ltkMpmHbMi6f3RJLLVgqF+vJSB0nRQGi2kqkC5itB1sExwPWS9EQh2u41IJPB7M7i5BGrTRi3W8bNJiufmSc06GIUWTs6i3a+T21lCNFq4h6ZQUkn8eh2OosV8LLqPZ7VE1/M8Pvaxj/H2t7+dzjr53vv6+vjkJz/J7/zO75DNHnMs/Lp8VlKpFCMjI0c87jgOBw8ePOLx1eRERHephep53qMs3VCMfd8nmUxG4qppGq1Wi1QqFQm0oii4rotpmti2TafTQVXVSIx934/+nUql6HQ6kUUcamAikUDXg/wA27ZpNpt4nhddIzwvFF2vm3dwrPdkZURXUYO0XdPAG8jjZQ3sjM7CeRq+DqlpycADQZytb6h4pop1YBGkpLOpn9aATmqqjX6gm16uqjgb+9FKgYtBum4QnZBIIC0jEFJVhUIJpMQ/dQzPVDEOLuKXK+A4QWywriEdF+ksnyYfi+7jWQ3R9TyPT3ziE9xwww3RY+F6oSgKl19+OZ/97GePVXjX7LMihGBgYIC3vvWtnHXWWbzkJS854rmVSoUvfOELj3r0nZ2d5cYbb0RKSafToXkCsezLcSKiq6rqo6xawzCi10OrFMB1XaSUCCHI5/MMDg6SSCSwLItcLkc+n+f+++8nm82SyWTYuXMnlUqFjRs3snHjRvbu3cvExASu60Ziq2ka7XY7un5oveZyOdLpNK7rUigUAKLPa+hq8Lp7R6Hl+1hrd9VFV8lkEOPD2EMZPFNBSKhu0KlsBaMiGL+jjvAk7UELJ6mg2pLMT6eQ2RRuLkFzzCJ3z2FkrQ6KQKRSyKSFn0+hFmr42SRuxsSYWMDvy6KUash2J7CGa3Xwuh806SMsCxIWslTGb7QQqvKkvt1YdB/PSouulJJ//ud/5q1vfeu6WbiPRQjB5Zdfzpvf/GZe+MIXHvXbjmH84/qshGL7hje8gbe85S309/fz6FIDR4fnedS6LrZ7772X733ve9x6663Mzs5GorESHK/oKopCKpXC931c132UxSulRNM0crlcZHVu376dc845h/PPP59UKoWu60xMTNBoNIDgi8f3fUzT5MCBA7RaLc477zxyuRyTk5Pce++9lMtl+vr6ADhw4EBkPYfXtywLy7IiH3D4JdZsNtF1nVarFfmfQ7G1bTv6Qujej1UWXSHQRkfwBntojacwSjZatY3Tm6S43cQ3BKPfm0eaGlJVqTwtg1nySO1ZoL25D2vfPHg+fqkc+IFNE5FKIjs2dDqIbIbi8zfgWoK+n1URngdCIGYLCF1HZlPYg2nqYwZSFRh1H73uYc7WYd8hxNgwlCp4i4UjLiEW3cez0qJ700038ba3vW3VrK0TIZvN8ulPf5orrrjiaE5f9c/Kn/zJn/A3f/M3jI+vfE3xyclJDh06xIc+9CEeeughdu/efcJjHq/oCiEiy1bTNHRdp9PpRL7ZRCLB4uIiuq4zPj7O8573PIaGhjBNM/LZho/5s7OzOI5Dp9OhUqmQz+dpNBqoqkoymSSZTLJv377I8rUsi6997WuMjY2hKAozMzPs378f2w6eipVuzkAo4qEot9vtR1nd4Tme50WujCcT3ROOXhCajkxaoAraPSqFM5N4VhKtAWobhAtubwonq9MY0RE++Kag/IwhyltVsoNj9N47G7gP0ikQCv5iAXQdYVk4pwzQGFFQbHB6LaxfzWFvGsCoJmif2s+hFxh4wzbmhMLQTxysuRai7SJsBxIWom0jRVw2eL2QUvLpT3+aG2644SkpuADVapXXve51KIrCy172snWbx+joKJ/4xCf43d/9XdLp9JO/4TgYHx9nfHycZz/72Rw+fJhbb72Vf/mXf1kR8T1WpJRRhICmaZGwJZNJdF2nXg9ckpdeeinPe97zogiC/fv3R5tumqZRKBSwLAtFUWi1WqTTaVRVxTAMMpkMk5OTZLNZNmzYwOjoKEII9u3bRyaTIZPJYNs2F154IU972tPYs2cPjUaDycnJyIL1PI+enh46nU5QO6br31264adpGkKIaA3LcWKWrhAIwwhqJXQnqI6NULxoGE8XqLbEMwV2RiA8sMo+VsGlPqZT3yBojzsgYeO/CZKHaoiWjajWkR0bOTbE4jN7KDzDR+kIEjMKA7+wSeyaxevPgSoonJOhvlFgLcLQf1VQHjmM9LxfJ0d0Q8iUvl68xcIRkyZiS/fxrISl6/s+n/vc53jzm98cPQI+lent7eUzn/kML33pS5c7bVU+KyMjI9x6660897nPPdq3rBjT09PccsstvPvd7z6u39OJWLq6rqNpGqZpRhto+Xyeer3O8PAw27Zt44wzzkAIQSKRYH5+nunp6cj/2m63yeVyKIpCtVrFtu1IxDVNo7e3l4mJCRKJBC94wQt45JFHkFJSr9cpl8u4rkulUuG0006jt7cXgHQ6zeTkJHfddRcHDhxA13UURYmiHjqdTuSe6XQ60Yaf67q4rhtutB3xnpyQCRhU+zLwW22EZaIOD9Lc2k87r5Cac1FtiZMOrq14kszBNsmdM/TsamJUAc1n0+Z5lLfM8cg7dGYvG6T+zFOYes3T2PNXSTa/7ldsPH0OL+8ifPAMBWkZqJUGbtpAqpCYk1hFH3sggchmgkpmvgyy1cLdz+rRhY7FrBy+7/OFL3yBN73pTSeF4AIUi0Ve85rX8N3vPrY41eoyPDy8boILgYV9ww03cPPNN3PZZZet2XWFEKhqENsf+lDDqIXBwUGuvfZaTj/9dCzLYmpqimKxiGEYWJZFo9Fg3759VCpB/XXbtlFVld7eXkzTJJ1OY5omxWIR13XZtGkTi4uLlEolJicnmZmZify3mUyGer2ObdssLi5y6NAhdF3n5S9/OX19fZEFW61WgSDKQQhBq9WKLN7QAg7XsxzH714QAiWXhd48JE2aY2kUV9LJa/gmFM7UMSoSNwlKB7QmaKUmeB76TAlrMUnVVji7Z5oLM/vp21znQ4M7KLaS7Bh+hEuyu3GkxvfLZ3Joppf8fpfUgQoYOr4qcBMqSMgc9kjMtaLQMNENnhamichnod0BTaP2jFESt/0U/JMrYy2fzwNQLpfXdR7Hyhe/+EWuv/76p6xL4UiUSiVe9apXccstt7Bjx45Vv97g4CC33HLLugnuUi6//HIuuugirrrqKu66666jSgI4EcKn7DAiIEyCqFarvOhFL4rErlgsUiqV6HQ6DA4ORuLqOA6pVIpkMsn8/HwUFpbJZNA0jYmJCYAoOWLXrl2Uy2XGxsYolUoAGIZBvV6nXq9TqVSwbZtEIkG73WZsbIwrrriCb3zjG1FcLwQGRTqdxnGcyJ0QriE8vhwn5NMVmTReLkl1Sxq1I1HbLlbJxagLGkMaVtnH1xSMWmDlIiV+qYyUkr77DJpDgzTONXGkioLPa8fuZk97hAtS+7GEw2Enza7SMIk9FumH52ChiH/qKG7WpJNXaYwKmkMaRjWDUZMkhxKoLR83pdIYUunZ20afKiMnZ0j/R5WTS24DNm3ahJTypBLdXbt28cY3vvGksXAfS6FQ4Oqrr+bmm29edeH94Ac/yCWXLNvGa00ZGBjgK1/5Ct/73vd4wxveEEU/rAZhFACAZQW9IB3H4fzzz2dwcJC9e/cyNjbGzp07WVhYQFEUyuUy4+PjpFIpcrkc7XabUqlEsViMQsBC//DCwgLj4+M0m0127dpFu90mk8mgKAqWZUU+Y9d1yeVyNJtNLMvCMAwqlQrtdpt0Os2OHTvIZDJ873vf48CBA6RSKQzDQFGUKK433EgL17Ucxy26QlXx00k8S6MxrNAakiAMPEMiDRkUHrcDX27mgEAvNoNsspC5Bfp/2cNd528hv61JJtPmoN2PjyCjtNlvD3Lz1LNYuGeEsbvbUCghXRfheHimiuJKZNeSN2oSo+ajtnwUz8co+yQmG6gzi0jPx7cd1FQSNZfF637DnSw8+OCD6z2FY6JSqfDud797Vf+zrgWLi4tce+21fOELX1i1R+4rrrjiyfzH60JPTw+vfOUrkVJy/fXXR4/VK034OK5pGslkkmq1imma9PX1RY/+YYJCMpmMNrASiQSDg4PRJtzk5GRk5RaLRRYXF3FdF8uyojCwgwcPomkaY2NBq77Qqja7tbtTqRTNZpOZmRk2bdpEqVTioYce4pRTTuE5z3kOGzdu5Je//CWTk5O0Wq3IX9xsNqNICs/zjurp4LhFV/oSZa6AUTYYmzXwM0ncnImbUKlt0GiOCjoDHlLz8Ewd0WjhFYooyWRQBMfzSf5qgex/jvAt7Uy2njXHwXYfHU9jrzHMz+qnMHHPBjbc2cbYNYlsNAP/sa4iVRASzCIM/KKNfv/eINtNKKD8+lvG17TA76xreMUSam/P8S435ij5zne+w5e//OX1nsaKMD8/z7XXXssPfvADzjrrrBUdO51O8573vCdyHz0Vufrqq2k0Glx33XWrMn6YGBHWU3Ach/7+frLZLNlslmQyydTUFGeccQa5XI56vU4ul2NhYYFKpcLi4mKUBrxp0yZqtRrT09P09PTg+z75fB7f97Esi3w+T6vVolwuRxlsYeRBJpOh3W7TarVot9vs37+fdrvN0NAQQ0NDVKtVvvrVr3L66aejaRo//OEPo4iFcBMtnD88uaV7fBtpYXxauYI3O4ecmYeH96FVOqi2T9/OFslpiVZVQAG9LpHtTrDB5bq/7m3Wsen7ZRt/f5rb5s7lJ3Mb2VMe5NbpC/j2w2cy8KCPuXsKuiLtN5sotTaKIzFLLtnDHuZEN/42FFw/6KcmhEDaDn6ng5JJx6191gDbtnn/+9+/3tNYUebm5rjxxhuj+M2V4oorruCcc85Z0TFXg5e+9KVcdNFFqza+4zgkk8mo2My2bdvIZrMsLCzQbrfp6+vj4osvZsuWLQwODlIsFpmamqJUKmFZFqOjoxw+fJg9e/YwMTERiWE6nWZ0dJSRkRGazSaJRIJEIkGr1WJxcRHDMMhms5GF6rou7XY7iooYHBwkm83S6XSYmJjA9322b9/OZZddxsDAAK1WCyDyIauqitmt5x3G+B6J47J0haYH7XIgENFkEiWVwDVV6mMGThJqm8HNeYimiuKCSCZQ0ylENhPE4RKItrHXo/ehzTzinYIUIBVQbcjNCrK/mPt15IH0kb5E1BoYxQTqYhULguI6jhP0TFNVRCYdVCArV6H7H0Wkkiiui9DUQHhPss20k4XbbruNhx56aL2nseJ84hOf4JJLLuEP/uAPVmQ8TdN4+9vfviJjrTbDw8PcdNNNnH/++Sueuh36RFutFs1mk1wux6WXXorjOJx++uls27aNBx98kJ///OfMzc1Rq9Wo1+v84he/IJPJkEgkyGazGIbBwYMHsW2bgYEBSqUShmEwMDCA7/v09PTQarUYHx+n3W5H7gRd1zEMI0rKSKVSLCwsAERREK7rRj87jkM+n+fss8/m7rvvjqIXQpfCY1OEj8Rxia50HZR8DmGZQXEaTcUdzlPdnMSo+9gZFS/tgeGTOKSRe6QZlGTs68HPJhH1epCiq2v4jSb5XzVQ7SRSESieRGv5JKYaUCh1O0p0w84SFtL1ELaLn05SPSOPZ4CvCjwriJDI76mjzZTA0HHP2YTScvFUBa0bxaD29eJ1b2zMymLb9oqmlz5V8DyP97///bz85S+PHidPhDDm9GRh+/btXH311Xz2s59d0XHD8LCl1uXw8DClUolEIsHs7CzT09OUy2UajQbVajWKVqjVapTLZQ4cOEC9XkfTtKicYy6XI5FIcPDgQdLpNPV6Hdd1o1Rex3GoVCrU63VSqVQU2xtGLaTTaWq1GpZl0W636enpiaxkXdc566yzuPvuux/lqsjlclSr1aP6/B//Rpqh4w3laZ07SnrnHOrDE/ROBJXGjO0jOGkTu0cht99H/+WBIGkhm8RLm2j9vdCYQmgaXq2GsvMRMpyG3WOCL7EOlpCzC/jdsmxhsXKRzSCrNZT5EjKdRK95LD5Xw9vURkrQJizye0A2mzSevYVOVqHvnhLYDn6zFbgYYmKOg+npaaanp9m4ceN6T2XN0TSNM88886isuGPB8zwcxyGdTkfXCMsvFovF6J632+0opbfdbrNp0yamp6eZn59/VDWxsNLYwsICw8PDLCwskMlkojFDa9UwDAqFAqqqsrCwgG3bWJYVbY75vk8ul6O3tzdKR9Z1nU2bNqHrOnNzc4yPj7N///5HVSELIxmWSziD4xBdYZqBxdlsoTQ6tPozLFwzRu7ACNn9TdRqG61mo7ZNeh+S9PznIXzbQRnsp74xg1TB0HswFotBUoWqBpltLQer2fWb6VrQeodu2x9dQ6RTOCN5dAhcB45D8r4Ko/I0phQLp9/FTUrmL8hgbN+Gakv67p7BOzwVZKYZBv5CEM0QE3OsTE9P85WvfIW/+Iu/WO+prAt/+qd/yt/93d+taFRKGN9ar9dJp9PRRlWz2aRQKLC4uEi1Wo0SIprNZhSGGCZANBoNdF2nWq0ihMBxHJrNZpSy2263SaVSUQQEBE9ktm2TTCZZXFyk0WgwOjoauRESiQSaptFoNPB9P7KIx8bGmJqaotVqcckllzAzM0O1Wo1C1ZaWkVyOYxZdadvIcHfu0DT9rof//CEWzxOUt6ZJTaeQAuw8JBYJxFNV8bNJPFOABKkIME0UwG+1ka6LWqoGPdV8H7obbhGehzPaS2PMIlc2EX46KO3ouCR3z5HYPo4zIPF7HOqbDPp/Jsk9XAbfR2gaytBAMIfp2W65x5XdFImJiTl2llblqlQq7N69G8dxGBgYiCqQhT5Y3/cpl8sMDw/j+z6JRIJUKoWUkrm5ORKJRCS0lmVF1cNs26ZYLEaREmEVs6XhXbquUyqVSKfTdDod2u021WqVU089lYWFBQzDiIqZV6tVOp0Op556KiMjI9HmW+gm6XQ6q7CRJhRQlEB8bQcxMcmQ7SD8cSpbJY0RgWqDl5B4hkB2hY/pBfSRNHZOozGqY87moVBGVuvBJlmtHqTuplPIRiPqdSY9D2wPJ2dgVjxw3KAnWsKEag1hWdh5UGoaelXhlG/V0OYrNM4YYu7lvfTuHiMxb2NMduNzVznL5reZJ/uwxQQs3Xz5bcb3fXRdj/yuP/3pT2k2m5x33nkA0aN/SJi4UCgUIt+u4ziRlVmtVqN03DAaotVqYds2rVYrSt8NNwQ9zyOdTkdiPDc3h2mabNiwAVVVo8gKx3GYnp4mkUjgui7Pec5zomuGro2wmlm9Xl8F94IiggiBWh0cF2XzBmi06P9JiUQxS2mrRmtI4vU5tAtGUPPWdWGgF63hIjWBa2n4CR1VCBTr18VyhKqC7QTugLCFRzKJ0FT0moO2WIeFAlI3wDSCnmlD/WQmJKlZD3OuhrpQxtnYTycfxPM6SUFuuoLoOEgIWgLFvdNWhZe+9KWceeaZ7Ny5c72n8pTGdV0++MEP8k//9E/rPZV1JfTDhqURZ2ZmmJyc5Cc/+QnPetaz8DyPYrEYFZoJC96Uy2VqtRr9/f1RFbBqtRp96Yf+WVVVo0iGcAPNMIzI6gUif2yj0SCVSkX+2zAeOJfL4boumzdvBoIaEeVyma9+9atRKJkQgmazGYWOrXj0gtA0yKaRowNIXaHZZ5Hcu4jSbJPe5ZB6RKcznMLOaYCHn0kganVEq4Pwkvjd6mNK047ifYVhBCFfltXtHmwHpR41LQhNM3TsnI7UMuiFUiCcjos0TZzeJP33LMB8AWGZtM8YozGi4yYEWgs6OYE7kEFpW6iuS/0ZG7G+ed+xLjvmKEin09FmQszyFIvF9Z7CuhJWFGu321F2WPjZ6XQ6HDx4kGw2Sz6fx3GcSIBDN0KlUsF1XTRNY8OGDei6jm3bLCwsRL7isF6vpmnU6/XIwu3p6Ylq74btekLrdGhoiE6nw+HDh+np6WHjxo3Mzc3RarXYuXMnU1NTfP/736dYLEbXDOccVk1zl7pGn4BjFl3fdtBqDdS2Da6LNq0GPczSJtojMwhfYu4/RGKgDz+fRnhdU9txUWsd7FMTSFUg2nZ3g0tHOm5QIUxKSCVQtSCUzG80Ee02pJMoth/4gl0X2WwidB3/lCFqG016F+tB+FpPFutAAeO+AkLT6Dz9VIQvUWttpKEhbTto6x4TE7OuKIoSPbqHtWjT6TRPe9rTGB4ejtro1Go1MplMJHCKopBOpxkZGaFWq5HL5SLrsre3N0rn9X0/SrDwPI9SqUQqlUJRFCqVCslkklQqRaPRiCzubDaLrutYloXneVQqFebn5xkZGYn8wJZlRdZt6BcO2/+Yphn5oJfjmEVXMXS8xQKi6x/BNFGqTaSuBp0e8jnodIIwr0Y7SP+t1RHNFoqhY9SzqB2JvbEX42B3HEvi5dPYfRa+oZDa5UKtjpKwAovXdlDbHkICmhZY25aJVARW2UMa3dfmC5BKBj8bOkahiTJXRDoOQtOQjSb64QLLfw/FxMSsNmHftvBR3DRNxsfH2b59O6qqUqlUotq1YQpv2EFiZmYmcjsYhsGGDRsolUoIITjrrLO477772Lt3L+l0Gtu2o5brYapuu92O6vem0+nIF1wqlejp6WF4eJjx8XFOO+00FEWJIhpe/OIX09PTw+7du/nqV78azSls9d5sNh9VxOdIHPPORxhyJZKJoK1Os4ksV1CLddANZKUGisCfmUPOLSIbzaDIuRNIneJItIaL0vGQqQQylcA5bYTFZ2SZvcjk4OWCheePIkYGuxt2TtAHTRXYOT3w43YC14Tda2GUbGR3Y09k0qAoYOiBWE8v4pXK+PUGslJFui7e7PyxLnldSaVSpFKp9Z5GDHDXXXet9xTWDcuyeNaznrVi44V+3LABZVhisVQqRYkKjUYDRVEif+74+DjVajUqkmNZFgMDA9i2HRXJeeCBBzhw4EDkggit0bAMY7gBt7CwQLFYpFqtMjo6GkVCTE1NUS6X8X0/cm+EKb+Dg4MIIbjooovYsWNHNLdw/qHFvPJVxmRXdJf67oQS1K11XYRpgBD49UZwjmkGFrCqIDWVVl9wycShFqLZRmoqi2cnMF82z+8NHiCrtfn66DmU20Pkvr2ItG3wfYxDRcRIPhBWTaWzqQ+pgjZfhVIliGgolUEoQYcIVQXHCWoxGDoimQw20Hz/pAoZu+KKK/A8j5tvvnm9p3JUnHvuuSddZbSj5Rvf+MaKjXX48OHIsjoZSKfTXHbZZdx+++0rMt6SXmKRqJbLZUzTJJvN0mg0ItGbn5/HsqwoXCtsWtnX1xdtlD388MPs2rUr6h5hGAbpdJpcLheFcyWTyagAjhCCYrHI/Pw89XqdkZGRKOqh0WgwNjaGpmns3buXdruNEIJsNhulKz/nOc/hrrvuiqIl0uk0xWIxyq5bjuMKGVMH+7E3DeAlNRK7BPg+frXbobfVQqSSCMNAtjtBrKyuISwTYTtoHYlZ7CCm5vA7HfxzttDeUeNVG37GNmuaWSfPlZsf5FO/9xyswumYP9qJ9HxkoYSma8h0ElFvYkxVMDQVFov4jeCbSEkFkRJ+u4PQu0tTBCJhIUf7URYryGYLTtI6rycDr3rVq/jMZz6z3tN4yvPjH/+Yffv28cxnPnO9p7IuhI/t4YZaaOnu2bOHsbExTj31VIrFIlJKms0mruvSarVwXfdRxW4mJibYuXMnv/zlL6MNuVwux9atW9m8eXOUAmwYBo1Gg8XFxaj1ztLQstBV4Xke2Ww2qkg2NjZGuVyOerCZpkmhUODAgQORVZ5MJmm1WlGPtpX36fbmkekklS0JPBOc9BiKLVHbHnrdQZ0p4RfLgQ/V8wILs9MBQ8cZ7yM11UbdN4VstVDyOfa9IsU1W3/IKcYizzQL7Fba5NUmr376j/nynks45cEUfqGI3xEoswvITaMw14CECTNFRC4LrTbStvGqHkJVUQw92KRT1SB1OJcG2a1wlk7BSVRTd6XbZa82YefVk61jxNHw2te+dkXHW1xcXNHxVpOwvOFKEvYTk1JGvdJ+/vOfk81meeUrX0mhUIiyysIkhLA9z1lnnUWj0WB6epqDBw9GSRFhZtvY2BjNZpNsNhs1ABgbG4vCusJ6C+Pj41E8bujy6O3txXVdhoeH6evrY3p6mmazyf79+xkbG2NmZoZ7770XIKrjG64ntHyX49gt3U4HFosM3N5GJkzam3qobNJJzQd90MRAsJvY2tKPcCW+oWAW2iiVJr6uYDw0gV8Nog1Kl2zm3Iv3oguPttRJCpWab7GvM4SCRL2whH3PKWg/qiBtG9/zUOaKIBTcoRxKTxql0QlShT0PVBUlncLfNILdY+ElVJyUQmLeQa8FO58ye3L5R3/wgx+s9xSOiYsvvpgdO3Zw2223rfdUVhRFUVbUpwnwvve9jx07dpwUYXZhxthKsTQbLRTepe3UZ2dnI3/q+Pg4vu9Tq9UisSwUCuzcuZOZmZmobc/WrVujJpXz8/NkMhlSqVTUCiisj5BMJnFdN3LteJ4XFboJxR2CFlm5XI58Ph8VOgqL6kxOTqLrOp7noev6o7qkPFmftGNPA3ZcBCArVagJEq6HNRkUkhG1Jn6xhDQMjFKGxXPSKC7oVQVF1zD3TONVqoG/9bQNzPyOzybF45HmAIN6lQdsi3vqW2l5Ogt2mpFMjZnzexh/KItstQPLOZUEw0VpOhSeniW5kCDl+Xj9aRaenqRyuofMuOAq4EBut0JuVxOpKUjfx983caxLjjkGhBC84x3v4Jvf/OZvVNbV1VdfzR//8R+v6Jg//vGPuf3223nxi1+8ouOuNOVymde//vVRe5uVIHQrhELreR6JRCKqndBoNCgWi5FPFgIx27p1K8lkkkajQTKZpL+/H1VVo3/39fVFPdESiQSlUom+vr4oRTf8O/T9+r7/qASJ4eFhbNum2Wxy4MABJicnGRoaIpvN0t/fz9zcHA888ACNRiOKA240GlGUROiyWI5jT45IJmGwN8gc0zXkzAKy0wn8pgCeh1+roeyfYqA1SG1bLsgkq9RwFwoIRaDkcxz8vTxbth6m7poUO0k+W76Iqf39aHmbF2zdjSokPoLmM1p0fr4Z475fIXwfezSP1BWMuTpqR9IYVilvGaR6hkO6vwyVBIm9JslZiVGXaB03CGcDZLN1xDbsMSvH+eefzzvf+U7e9773/UYIb09PD29605uigikrheu6vOMd72Dbtm1RxtNTkTvuuGPFn7jCughhQkGxWIwKzpRKJR5++OGoQhgE0ROpVIpt27YBQQPRkZERHMdhYWGBiYmJqONvtVolmUxi2zb9/f1MTU1FGWeqqqLrOrquI4RgYmIiavteq9WiTh5h14mw31q73WZmZoaHH36Y3bt3R4kZYaxxmBQRWr7LcRwZaSq0OjRPH6Bwpk5yrpfs/hb6/lnkkvRaoQdhW2bZhUotKCReb4CUONvHaJ7qMFdLkzJtCuU0iXtTbPuvGtXTUvz772/naWOz1G2DbKbJgZfl2aRuI7FvkcaYha9Dz4JCbm8DtdZm9tJ+rCmd9H9kGdvXQCmXELYDioIzkkdqCmqlhVzhIswxT4xhGLzrXe/i4MGDfO5zn1vv6ZwQ+XyeT37ykzz72c9elfEffPBBrrnmGu68886o88BTiW9+85u8/vWvX/EvzzB6AX5d4jEUxHa7TS6XwzAM9u/fz8aNGxkeHsZ1XQ4fPhxZqhAUWQ+TKVqtFq1Wi1qtRjabpVarReUcw2SHZDJJoVAgn8/T399Pq9WKEiUSiQT1ep2hoSF6enrIZDJ4nsf8/DzNZpNKpcKuXbui7LlMJoPrulHLICDqubbs2pc74QXKHz7uoJLJBK3XAfu0Qdp9BnZaQWtL8vdNg+shWy3sczZRGzfp/88p/GyS+pYc2XsP09k6zPTzLFrjDmraRc6bbPiBR2rnHP5CASWboXXmGBOX61hjdYSQ2B0d03JoHwhKQ2YeURi+u4LoeCi1BvbGfvTFepAc4ctuuFq3d9FgP1LXoFjG78bqPpZ/97+0fGDd0hsmxPJ39DcEKeVR35PwLY994f777+fKK6/k0KFDKzSrtSWbzfKpT32KK6+8cunLK/5ZsSyLf/iHf+C66657Un/gWvLd736XV73qVRQKhWXPO5bPSnhPwuLfQGQZptNpMpkMrVaLoaEhduzYgW3bqKqKZVlMTU1FG121Wg3f96NIhMXFxSi29tChQwwNDVEul1FVlYGBAarVKrt37yabzeJ5HmeddVYkjqHIh6Fquq5HvdHCbhEQWPwPPPBAuOYolti27SgNOLTcq9XqEe/Jsft0W63gW8/z0Iol0kIgRgZpbuvHz6VwsxZ6oUFz0ED4Ejo2bk8ftTGV5ktOoXSOj9RcUMD8ZYINP6ghHt6PTFiBUFom1s8O8LQ9SWZesoHSM1yE7uPpHr3bC1TqCZQ9KZRqi86GHnRdRd89iV+uID0PJRxHCPC8ILHCl4EvOW7Ts6ZccMEFfOUrX+EVr3gFhw8fXu/pHBPpdJpPfvKTjxXcVaHdbvO2t70NVVV53ete95So1nb77bfz6le/+kkF90TQdR0pZRTSFfphU6kUc3NzfPe732X79u0MDAxENW3Dfmeha0JVVcrlMtVqNSqMc8YZZ7B///4o9Evt9ky0bZtyuUwymYyuPTw8zNzcHIqi0Gw2qdVqqKrKaaedBgR+5C1btvCv//qv3HfffaiqGoltrVaL5r80KeLJ2hodV+cIpb8XhMCbmQuSDwplhNtHc0OG9MPz0GqTmk0hXAmGTmvAoLrVQ+mzUQDDdGgVEqSnJKLV7aBpWfiVKmgq/qljeKZKfp+NVdJpDOtUz1QCq7dsMjjrw0IB+6wB7HyajO2iuC709+D1pcGXaIcXgrKSuh64PWLBXRcuuOACbr31Vq688kqmp6fXezpHRSqV4qabbloTwQ3pdDq85S1v4eDBg1x//fWMj4+v2bWXUq/Xuemmm/j7v//7VQ1pC9upq6qK7/tRQkO1WiWbzeL7PjMzMxQKBc455xxOPfVUkskkhmFgmmbU7WFp/dpGo0E+n8c0TQYHBzFNMzpX0zSy2SzpdJrNmzeTSqXodDoMDAyg6zoPPfQQzWaTwcFBLMuiXq8zOztLuVzmM5/5DIcOHYoy4cIQtqUZdaH7JdwcXI7j65Gma8iUhVJJgvQRuQxOWkVv+MjFIqgq+nQFmTSRlonwgaxLNtPE8VSaB7JkDyu4lkSpN/EhKGTTasH0HGIajN4ekBKrWiMPDJ57GnPPytJfkiSn6jAyiGr7GGUHUa4FdXeFQHgSYbt4o30A+JaO9qup41lmzApx0UUXccstt/BHf/RHzM7Orvd0liWdTnPjjTdy1VVXPWk650rTbrd573vfyx133MGXvvQlRkdH1/T6tVqN6667jltvvfVJ/ZInSthlYek9DgvXhMkIYXuen/zkJ9i2zYUXXoimaRSLRfL5fFR8ptV9+g5jdQ3DiNwLoQUbuhlc12VsbCzq8tvX10cymYzies8//3yKxSJ79uzhnnvuYWZmJirIA0HXiaVhbuGG2tKmlCteZUx6HsL1cDImen8P7VN6KG0zSSz6JHfP4XcnpRQFipeFjo21kEGUTUp+BmvC4LR/byBVQeW0RFBfodMJhFdVg0I5yQStbUNY0zVwHPxCEWPfDOOzSTB0/IRO8Rl9qDYkDtbwq7WgItn+JiKRQEkFu8wylwnig0+iZIjfVJ73vOdx880388pXvpK5ubn1ns4Tkk6n+djHPsY111yzrvO4++67ueqqq/jgBz/IBRdcsCbuhttuu42//uu/XrMU7rAObmg9hhECuq5HBcTDDhC+7/Pggw9Sq9UYHR1lbGyM3t7eyC8cWruh4FUqFarVKv39/czPz1MoFBgcHGTTpk3Mz89HCRaTk5OYpsnFF1/Mhg0baDQa7N27l4mJCe68805arVYUFhYWxgnTkEORDQU2DH1bnTRgKZGNFmrDoX5GP5VNGpkpj9xPpoPYXUAxzSATTBFQq6OVslhzCYxfGYzcuYg8PINQFfJyM+3NfVi2E7RNL5ZgdIjGaXmEL4PaDJ6HunEcryeFb+m0Bk06OQUnLTAqErcvgVHLo2QzgQ+3062roKowX8BvNp9w8yxm7bn00ksj4Z2ff2oVHkqlUvzjP/4jr3nNa9Z7KgD86Ec/4vnPfz4f+MAHuOqqqxgcHFyV65RKJf7H//gffPnLX17zGr9hjCwQRRmElmNYmFzX9aiG7u7du9m3bx9jY2M885nPZHh4mI0bN0alHEOhrtfrzMzMkM/n6esLnnir1SobNmzAMAxGRkYAeNGLXsQDDzxAX18fk5OTfOtb32LPnj1RPLKiKFGZx3q9TqPRiLoOL808C8U//IJ4siekY45eCNGGh0DTcMf7aA1bZH4xjz+3EMTrjg1SPT2LakvSu4uIap3muRuwZht4ST0oUrNYDIrh9OUDt0C1gb+wiNLXi0xaQciX6+FXayj5HDgOMpehenYfZtFFr7Sx+xKB/7bhMndhCqlAcs7HKnpobQ/j4Um8uSf/zx1HLzyelYheOBJ33HEHV111FQsLC8d4idXBsiw++tGP8rrXve5oTl/Tz4oQgjPPPJM3vvGNvOY1ryGbzZ7okJGg3X777Xz84x/njjvuOKHxjid6YSmapj2qtf1S4Qp9uGH9Wtu2o1jYbDbLC17wAl74wheiaRo//elPqVQqDA4OMjU1haqqXHzxxUxMTDA1NUUqlcJxHM4991xarRY/+tGPuPPOO+nt7WVxcTEKLRNCoGnao3qqtbvhsOGGXOiLDoulhwkeoZ4ud0+OW3SFbqCODuEN5BBtF6VQBimxt4xQPj1Bp0dgFiW5Ax30n+5FyeeQtRqitwevJ4XYPRHURkgmgo0uVQ3Sg7uhXtK2o5q9or83uKjtBHUUpmYRqVRQaGeoF2G7zFzaj29Az68crJkmSsfB3ztxVBXFYtF9PKspuhCkN19zzTUsLi6uuv9wORKJRBSudZSP8evyWdE0je3bt3PDDTdw7bXXRqJwNPi+H23u3HrrrTz00EPceOONUdfcE+VERReIUnSXWomu62KaZmT9hqKpqmpUYMYwDMbGxjj77LOjBJNqtcqBAwdYWFjgvPPO4+DBg1QqFTKZDLt370bTNDqdTtTWPWzjHt6L0J0AgZ87THoI5xZGKoSvh+K79HO8KqKLoqL25ILoBdPAH+pl/sI8pXN8tJqCXgs6/yYWJAP3VxCHZ4N2PM0WIpcBzwdVASnxi+WgdoIQKJk0Mp2EUgXZCsQYz0Nk0sjhPkTLhvnFQHQVBRSBO9JDeWsKre2T+0UBFot4pcpRRyzEovt4Vlt0IShv+K53vWvdqpJdcskl3HjjjWzZsuVY4mPX9bMSFmm57LLLjjp9+Nvf/nZUC3hqampFhHYpKyG63WNRKFe4KRUKXehbBaJSjUD0WijQ4WN+mIq71AJd2h49HCOM9Q3TksP+amHSw9JuxEu7AofFbcKxH6ujqyO6XbThIfz+HlobM9jZ4IOruJLqRhXFBeGBWfHpvW8B5heR7Q7SdZF+0IhS6BoimQha9SgqIp0Ex8Xvy6LMl/DrDZR0Cm+4D7vPwlxsoVQa0O4E7doB+4xxKptNtJak5/55/InJY6qZG4vu41kL0YXAknjTm97E5z//+eN5+3GRyWTYsWMHH//4xxkeHj7Wt8eflcewUqLbPQ782uUQblAB0SYbEAlgeF6n08GyLKSUj4qTDf3Guq5jmma0Eea6bpQBF4Z9hZ0gQqF2HCe6fmhVh9cK6zYsuQdHfU9OWHTV/r6gDoOhI1JJWtuGaAzrJOddrJk6ou3gDGZQWw7qbAm6jnN76yilbRZShb5ftlB/8UjgSshlg4piGwZxcibGQgO7P4XaclFbDkqhijSDDBbhesh6A3rzyKQJUiIOzeIdY7RCLLqPZ61EF4LY0K997Wu8/e1vX/WQslNOOYW//du/5dWvfvXxDhF/Vh7DSopu95zIn9odP/o7tIDDMK1QDMMuv6ZpRskQQFRAJ5lMRiFouq5HYhqWjlxaZhJ+7ZIJ5xJGKoQbdo9Z/zHdk+OK012Kt1hASSYRG0ZonJLFTSr03V8I4m0TQTseY98M7imDtLYPo7Y9jKkSbkqjMSrI7/PRZ4L2GBIQUiKH+7DzJr6h4OYSeJaC+YtJZCf4dhHpFCKZwE8nEaaBqNQQrQ5oKm63dmbMyUM6nebVr341W7du5cMf/jC33HLLqlznuc99Ll//+tejHe2YpybhY/vSWrWh2yEUwfCRPrRkw+SKTqcTxdWGm25A5L8Fokw2IOoKEVq+SzPMQrdDuIF3BDfCMa/vhC3d7ipQ+3qDwuXtNn6rHbRVTybwxwdQam1kwqB0dg7XFPT/vIZSaVI7awCt6QXxvQsFpOOinDJGY/tANLTecDH3zuGXysGlLBPZ7qD05HE29qPWO4iWjWjbeDOzxxUeFlu6j+c4LN2YmJij4IQtXQCkxFssBKFfhgG+RLZaSNtG6XSCTa9SFW1LFqkIlEYHP5sgNdlELdTAdlCyGfxyBf/wNImDUyipBJhm0PyyXg822LpfECKtBe3fK22YC1IV3cXVyxGPiYmJWSlWRnRDpATPQ0kmQRGgG4hUAjo2cmwAs+yidlREuYZaC8z7zpYhfE1gPTgRWLGtNkg/6B7s+UE7dVUFTYN2G2EFbdml7yMnJvGbzeC6MTExMScBKyu6EEQmhDG2ho4z2hNU/PIlermNtXMRv1ZHJJOIVAKpCIyKHYSSGTpC01B68oFV2+lEDS9lqRL0P7OdwIWgqrHgxsTEnHSsuOgC+O02tNsorou2px1EJBRLSM/Hd5ygOeRAD3KhhLVnBgCZSjFz1VZG7igg54uQz+AOjKAfLuDNLyBQEYYepPX6Mq4aFhMTc1KyKqIbIlstZKuFmsuCUPAb9SA2F/CTBmJsAKYWkLUaytAAel2C4waZa66LVizjVetINwjRkJ4Xi21MTMxJzcpELxzVlbqb4UJBscyg4Hg+F2SiQVCEfJ0K08TRC48njl6IiVkdVtXSfRShuEsv8MXCURWjiYmJiflNYv37gsTExMT8FhGLbkxMTMwaEotuTExMzBoSi25MTEzMGhKLbkxMTMwaEotuTExMzBoSi25MTEzMGhKLbkxMTMwasmxGWkxMTEzMyhJbujExMTFrSCy6MTExMWtILLoxMTExa0gsujExMTFrSCy6MTExMWtILLoxMTExa8j/D1366ANatlPHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Computing worst-case/best-case SNR using the corpus callosum...')\n",
    "from dipy.segment.mask import segment_from_cfa\n",
    "from dipy.segment.mask import bounding_box\n",
    "\n",
    "threshold = (0.6, 1, 0, 0.1, 0, 0.1)\n",
    "\n",
    "CC_box = np.zeros_like(data[..., 0])\n",
    "\n",
    "mins, maxs = bounding_box(mask)\n",
    "mins = np.array(mins)\n",
    "maxs = np.array(maxs)\n",
    "diff = (maxs - mins) // 4\n",
    "bounds_min = mins + diff\n",
    "bounds_max = maxs - diff\n",
    "\n",
    "CC_box[bounds_min[0]:bounds_max[0],\n",
    "       bounds_min[1]:bounds_max[1],\n",
    "       bounds_min[2]:bounds_max[2]] = 1\n",
    "\n",
    "mask_cc_part, cfa = segment_from_cfa(tensorfit, CC_box, threshold,\n",
    "                                     return_cfa=True)\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "region = 40\n",
    "\n",
    "# visualize the masks and ROI\n",
    "fig = plt.figure('Corpus callosum segmentation')\n",
    "plt.subplot(1, 5, 1)\n",
    "plt.title(\"CC\")\n",
    "plt.axis('off')\n",
    "red = cfa[..., 0]\n",
    "plt.imshow(np.rot90(red[...,region]))\n",
    "\n",
    "plt.subplot(1, 5, 2)\n",
    "plt.title(\"CC mask\")\n",
    "plt.axis('off')\n",
    "plt.imshow(np.rot90(mask_cc_part[...,region]), cmap='gray')\n",
    "\n",
    "from scipy.ndimage.morphology import binary_dilation\n",
    "mask_noise = binary_dilation(mask, iterations=10)\n",
    "\n",
    "mask_noise = ~mask_noise\n",
    "print(mask_noise.shape)\n",
    "\n",
    "plt.subplot(1, 5, 3)\n",
    "plt.title(\"noise mask\")\n",
    "plt.axis('off')\n",
    "plt.imshow(mask_noise[:,:,40], cmap='gray')\n",
    "\n",
    "plt.subplot(1, 5, 4)\n",
    "plt.title(\"mask\")\n",
    "plt.axis('off')\n",
    "plt.imshow(mask[:,:,40], cmap='gray')\n",
    "\n",
    "\n",
    "plt.subplot(1, 5, 5)\n",
    "plt.title(\"data\")\n",
    "plt.axis('off')\n",
    "plt.imshow(data[:,:,40,40], cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize raw data for metric computation\n",
    "data_normalized = data - np.min(data, axis=(0,1), keepdims=True)\n",
    "data_normalized = (data_normalized.astype(np.float32) / np.max(data_normalized, axis=(0,1), keepdims=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_signal = np.mean(data_normalized[mask_cc_part], axis=0)\n",
    "noise_std = np.std(data_normalized[mask_noise, :], axis=0)\n",
    "mean_bg = np.mean(data_normalized[mask_noise, :], axis=0)\n",
    "\n",
    "mean_signal_denoised = {}\n",
    "denoised_noise_std = {}\n",
    "denoised_mean_bg = {}\n",
    "\n",
    "# normalize denoised data for metric calculation\n",
    "for k, v in denoised_data.items():\n",
    "    mean_signal_denoised[k] = np.mean(v[mask_cc_part], axis=0)\n",
    "    denoised_noise_std[k] = np.std(v[mask_noise, :], axis=0)\n",
    "    denoised_mean_bg[k] = np.mean(v[mask_noise, :], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "raw [SNR] mean: 5.1141 std: 2.4988\n",
      "raw [CNR] mean: 4.6567 std: 2.4976\n",
      "=================\n",
      "n2s [SNR delta] mean: -0.5301 std: 0.7778 best: 0.5438 worst: -2.9176\n",
      "n2s [CNR delta] mean: -0.5535 std: 0.7710 best: 0.5077 worst: -2.9232\n",
      "=================\n",
      "n2n [SNR delta] mean: 0.8060 std: 2.0059 best: 4.0267 worst: -3.2664\n",
      "n2n [CNR delta] mean: 0.7942 std: 2.0026 best: 4.0074 worst: -3.2668\n",
      "=================\n",
      "dip [SNR delta] mean: -3.0587 std: 1.6048 best: -0.2501 worst: -6.3830\n",
      "dip [CNR delta] mean: -3.0380 std: 1.6036 best: -0.2839 worst: -6.3900\n",
      "=================\n",
      "p2s [SNR delta] mean: 0.1042 std: 0.3641 best: 0.8540 worst: -0.8332\n",
      "p2s [CNR delta] mean: 0.0748 std: 0.3623 best: 0.8237 worst: -0.8545\n",
      "=================\n",
      "our [SNR delta] mean: 1.8284 std: 1.6969 best: 4.9025 worst: -1.7451\n",
      "our [CNR delta] mean: 1.7486 std: 1.6949 best: 4.8205 worst: -1.8113\n"
     ]
    }
   ],
   "source": [
    "# metric calculation starts here\n",
    "\n",
    "SNRs = {}\n",
    "CNRs = {}\n",
    "\n",
    "SNR = mean_signal / noise_std\n",
    "CNR = (mean_signal - mean_bg) / noise_std\n",
    "\n",
    "SNR = SNR[sel_b]\n",
    "CNR = CNR[sel_b]\n",
    "\n",
    "SNRs['raw'] = SNR\n",
    "CNRs['raw'] = CNR\n",
    "\n",
    "print('raw', '[SNR] mean: %.4f std: %.4f' % (np.mean(SNR), np.std(SNR)))\n",
    "print('raw', '[CNR] mean: %.4f std: %.4f' % (np.mean(CNR), np.std(CNR)))\n",
    "\n",
    "for k in denoised_data.keys():\n",
    "    SNR = mean_signal_denoised[k] / (denoised_noise_std[k]+1e-7)\n",
    "    CNR = (mean_signal_denoised[k] - denoised_mean_bg[k]) / (denoised_noise_std[k]+1e-7)\n",
    "\n",
    "    SNR = SNR[sel_b]\n",
    "    CNR = CNR[sel_b]\n",
    "    SNR -= SNRs['raw']\n",
    "    CNR -= CNRs['raw']\n",
    "    SNRs[k] = SNR\n",
    "    CNRs[k] = CNR\n",
    "    print(\"=================\")\n",
    "    print(k, '[SNR delta] mean: %.4f std: %.4f best: %.4f worst: %.4f' % (np.mean(SNR), np.std(SNR), np.max(SNR), np.min(SNR)))\n",
    "    print(k, '[CNR delta] mean: %.4f std: %.4f best: %.4f worst: %.4f' % (np.mean(CNR), np.std(CNR), np.max(CNR), np.min(CNR)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz8AAAE/CAYAAABoyn1qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr50lEQVR4nO3deXwV5b3H8W8gJITFBJBgpZfIIigE8GUoi9SkFNuyBxBtUMmFKmUvVhQiBFoMJIBgKuYlXL0oF2tdEAStQRHUpFKKfeWWi4GyqJAGkAbIYhLgZJv7h80pEbJwksycw/N5/xXmnJnnNw+T8+Q788wcP8uyLAEAAADAda6J0wUAAAAAgB0IPwAAAACMQPgBAAAAYATCDwAAAAAj+Ez4sSxLLpdLPJ8BAAAAQHVqyg3+DtTjEZfLpYMHD6p79+4KCAhwuhwAAAAAXqikpERHjx5Vr1691Lx58yqv+Uz4KS0tlSQdPXrU4UoAAAAAeLvS0lLfDT/NmjWTJK78AAAAAKhW5ZWfyvxwOZ8JP35+fpKkgIAABQYGOlwNAAAAAG9WmR8u5zMPPAAAAACA+iD8AAAAADAC4QcAAACAEQg/AAAAAIxA+AEAAABgBMIPAAAAACMQfgAAAAAYgfADAAAAwAg+8yWnAADPpKenKy0tzeP18/PzJUkhISEerR8VFaXIyEiP2wcAoKFw5QcAUKOCggIVFBQ4XQYAAPXGlR8AuM5FRkbW68pLQkKCJGnx4sUNVRIAwGH1mRVQ3xkBknOzArjyAwAAAKDOfHlGAFd+AAAAAMPUZ1aAL88I4MoPAAAAACMQfgAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABiB8AMAAADACIQfAAAAAEZwNPysXLlScXFxTpYAAAAAwBCOhZ+9e/fq7bffdqp5AAAAAIbxd6LR/Px8JScna/r06Tp8+LATJaCBpKenKy0tzeP18/PzJUkhISEerR8VFaXIyEiP2wcANCynxwWJsQFA9Ry58rNkyRL9+te/1g033OBE8/AiBQUFKigocLoMAICXYFwA0Jhsv/KzefNmfe9739OgQYO0devWa14/MzOzEaqCp1q2bKkRI0Z4vP4bb7whSfXaRkZGhsfrAqhdYWGhJH7XUDfeMC5IHK9AY/LlccH28JOamqqzZ88qOjpaBQUFunDhghITE7Vw4cI6rR8eHq7AwMBGrhJ2SU1NlSRFREQ4XAmA6vB7CjtxvAHez9t/T10uV7UXTGwPPy+//LL7561bt+qzzz6rc/ABAAAAAE/xPT8AAAAAjODI094qjR8/XuPHj3eyBAAAAACG4MoPAAAAACMQfgAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABiB8AMAAADACIQfAAAAAEbwd7oAAEDtNm3apKysLEfarmw3ISHBkfbDwsIUGxvrSNsAgOsL4QcAfEBWVpaOHz2sjq2a2d52K6tcklRy+kvb2z5VVGp7mwCA6xfhBwB8RMdWzTSnT6jTZdjquQM5TpcAAF7J5BkBkuezAgg/AAAAgI8xdUaAVL9ZAYQfAAAAwAeZOCNAqt+sAJ72BgAAAMAIhB8AAAAARiD8AAAAADAC4QcAAACAEQg/AAAAAIxA+AEAAABgBMIPAAAAACMQfgAAAAAYgS85NdymTZuUlZXlWPuVbSckJDjSflhYmGJjYx1pGwC8lZNjA+MCgMZE+DFcVlaWDn9xTM3atnak/fJ/HYFf5p6xve3S3ELb2wQAX+Dk2MC4AKAxEX6gZm1bq91PBzhdhu3O79zndAkA4LVMHBsYF4DrH/f8AAAAADAC4QcAAACAEQg/AAAAAIxA+AEAAABgBMIPAAAAACMQfgAAAAAYgfADAAAAwAiOfM9PSkqKduzYIUmKiorS/PnznSgDAAAAgEFsv/Lz5z//WZ9++qnefvttbdu2TQcPHtSHH35odxkAAAAADGP7lZ/27dsrLi5OAQEBkqSuXbvq9OnTdpcBAD4lPz9feUWleu5AjtOl2OpUUana5Oc7XQYA4Dphe/i59dZb3T+fOHFCqampev311+u8fmZmZmOUZazCwkKnS3BUYWGhMjIynC4DqNWlS5ecLsExly5d4vfUZiaPDYwL8BWFhYUKdLoIB3n6u+rIPT+SdOzYMU2bNk0LFizQLbfcUuf1wsPDFRho8n91w0pNTVVObrHTZTimdevWioiIcLoMoFapqakqqSjWnD6hTpdiq+cO5Cjgppv4PbWZyWMD4wJ8xR/+8AedMnBGgPSvWQEtyqv9XXW5XNVeMHHkaW8ZGRmaPHmy5s2bp3HjxjlRAgAAAADD2H7l5+uvv9asWbOUnJysQYMG2d08AAAA4PNCQkLU4sJ542YESP+aFRAS4tG6toefDRs2yOVyacWKFe5lMTExmjhxot2lAAAAADCI7eEnPj5e8fHxdjcLAAAAwHCO3PMDAAAAAHYj/AAAAAAwAuEHAAAAgBEIPwAAAACMQPgBAAAAYATCDwAAAAAjEH4AAAAAGIHwAwAAAMAIhB8AAAAARiD8AAAAADAC4QcAAACAEQg/AAAAAIzg73QBcFZ+fr5Kcwt1fuc+p0uxXWluofKbNHe6DADwOqaODYwLwPWPKz8AAAAAjMCVH8OFhITofMUltfvpAKdLsd35nfsUEhLidBkA4HVMHRsYF4DrH1d+AAAAABiB8AMAAADACIQfAAAAAEbgnh/AAenp6UpLS/N4/fz8fEnyeG56VFSUIiMjPW4fAADAFxF+AB9UUFAgyfPwAwDwPvU5MVbfk2ISJ8ZgBsIP4IDIyMh6DTAJCQmSpMWLFzdUSQAAH8ZJMaBuCD8AAABeoD4nxjgpBtQNDzwAAAAAYATCDwAAAAAjEH4AAAAAGIF7fgD4DNMfEX6qqFTPHcixvd3CknJJUuuApra3faqoVJ1tbxUAcL0i/AAwhi8/DSksLMyxtouysiRJ7W62v4bOcnbfAXg/HhGOa0H4AeAzTH5EeGxsrGNt+3K/AUBNfPmkmGTmjACpfrMCCD8AAADwWaY+ItzUGQFS/WYFEH4AD2zatElZ//rFd0Jl25Uf2nYLCwtz9EoEAHgjJ8cGxgXzMCPAM4QfwANZWVk6euxLtbzhRkfar/ALkCSd+meB7W0Xf3PO9jYBwBc4OTYwLgB1Q/gBPNTyhhvVZ+BYp8uw3YG/bHO6BADwWiaODYwL8CV8zw8AAAAAIxB+AAAAABiBaW+AB/Lz81X8zXkjL/UXf3NO+YGW02UAAABcM0fCz7vvvqt169aptLRUkydP1oMPPuhEGQBsxlPyeBoScD0z9cRYfU+KmfyUPImxwW62h59//vOfSk5O1tatWxUQEKCYmBgNGDBA3bp1s7sUwGMhISEqdvkZd1Or9O2NrSEhwR6tm5WVpeNHv9RNrdo1cFV108JqJkm6eDrf9rbPFJ23vU0A8AVOjg1OjgsSY4MTbA8/f/7znzVw4ED3N+n+7Gc/0/vvv6/Zs2fbXQoAB9zUqp1+0WeM02XY7qUD7zhdAoBGZuqJsfqcFKvE2AC72B5+cnJy1L59e/e/Q0NDdeDAgTqvn5mZ2RhlGauwsNDpEhxVWFiojIwMj9YzWX36zeQbDT3tN6dVHu++WDs8Y/JnXH1+T+k3z/uNscG3Pl99eVyw/VizrCvnhPr5+dV5/fDwcAUGBjZkSUZLTU1VTm6x02U4pnXr1oqIiLjm9VJTU/XNBfu/SM5b1KffLhbmN3xBPsLTfnNaamqqJPlk7fCMyWNDfX5PTR4b6ttvjA2+9fnq7eOCy+Wq9oKJ7Y+67tChg86d+/c3Aefk5Cg0NNTuMgAAAAAYxvYrP3fddZeee+455ebmKigoSDt37nT0CRuAp4q/OefYE31KXBckSQGBLWxvu/ibc1KH+s3tBgAAcEKdws+lS5eUlZWl7t27y+VyqXnz5h432KFDB/36179WbGysSktLNWHCBPXp08fj7QFOCAsLc7T9rKx8SVLHDt+zv/EOwY7vPwB4K6dOjHFSDKibWsPP/v37NXv2bPn7++v1119XdHS01q1bpzvvvNPjRkePHq3Ro0d7vH5jSE9PV1pamkfr5ufnS5L7CXaeiIqKUmRkpMfrw15OP4+/8mrp4sWLHa0DuJ7VZ1yQGBtM5OSJIU6KAXVTa/hZtWqVNm7cqMcff1w33XSTVq1apeXLl2vLli121OcTCgq+vbmxPgMcAOD6wthgHidPjHFSDKibWsPPpUuXqnwBaVRUlJKTkxu1KCdERkZ6fHaNDxygbvLz85VbdN7I7zX4uui82uY7XQWuRX3GBYmxAQC8Ua3hx9/fXwUFBe7HUX/11VeNXhQAAADMwIkxp6swS63hZ8aMGXrooYd07tw5PfbYY9qzZ4+eeuopO2oDcJ0JCQlR4AUZ+y3eQUx/AgDAUbWGnyFDhqhLly7as2ePKioqNHPmzCrT4AAAAABPcWIsxOkyjFJr+MnPz1dwcLBGjBhRZRk3cAIAAADwJbWGn4EDB7rv96nUvn17paenN1pRAAAAANDQag0/hw8fdv9cWlqqnTt3VlkGAAAAAL6gybW8uVmzZho5cqT27NnTWPUAAAAAQKOo0z0/lSzLUmZmpr755pvGrAkAAAAAGlyd7/mxLEuS1K5dOy1atKjRCwMAAACAhnRN9/wAaBjp6elKS0vzeP2srCxJ//4G+WsVFRVVr2+uBwAA8EXVhp+XX365xhWnTJnS4MUAqJvg4GCnSwAANLD6nBir70kxiRNjMEO14efo0aN21gEYJTIykgEGANBgOCkG1E214ScpKcnOOgAAAIzGiTGg8dV6z8/f/vY3vfDCC7pw4YIsy1JFRYVOnjypTz75xIbyYIfS3EKd37nPkbbLL7okSU2DAm1vuzS3UGp7k+3tAoAvcGpsYFwA0JhqDT/x8fGKjo7WBx98oJiYGO3evVs//elP7agNNggLC3O0/co5ymFODDZtb3J8/wHAGzn52ci4AKAx1Rp+/Pz89Mtf/lJ5eXnq0qWLxowZo4kTJ9pRG2wQGxvraPuVN2YuXrzY0ToAAP/m5NjAuACgMdUaflq0aCFJ6tSpk44dO6aIiAiVl5c3emEArk9nis7rpQPvONJ2UckFSVKrgBa2t32m6Lw6K8T2dgHAFzg1Njg5LkiMDU6oNfz07dtXjz76qObOnatp06bpxIkTatq0qR21AbjOOD2dJCerQJLU/uYQ29vurBDH9x8AvJGTn41OjgsSY4MTag0/OTk5uu2229S5c2ctWrRIe/bs0erVq+2oDcB1hmmWAIDvYpol7NSktjcMGjRIH3/8sX7yk5/o8OHDmjp1qrp06WJHbQAAAADQYGoNPzExMXrzzTe1fv16FRQUKCYmRrNmzbKjNgAAAABoMLWGn0qXLl1SSUmJLMvinh8AAAAAPqfWe35eeuklvf322yopKdGECRP05ptv6sYbb7SjNgAAAABoMLWGn4MHDyo+Pl4DBgywox4AAAAAaBS1hp81a9bYUUe9bdq0yf2t0HarbLfyiSFOCAsLc/xJWgDgbUweGxgXAOBKtYYfX5GVlaUvvvhCbdq0sb3tynugzp8/b3vbkpSXl+dIuwDg7UwdGxgXAODqrpvwI0lt2rTRPffc43QZttu1a5fTJQCA1zJxbGBcAICrq/PT3gAAAADAlxF+AAAAABjhupr2BgC4Unp6utLS0jxev7437kdFRSkyMtLj9gEAaCiEHwBAjYKDg50uAQCABkH4AYDrXGRkJFdeAABV1GdWQEM8yt+pWQGEHwAAAAB15sszAmwPPxkZGUpMTFRZWZlCQkKUmJiojh072l0GAAAAYCxTZwXY/rS3J554QsuXL9f27ds1evRoLVu2zO4SAAAAABjI1vBTUlKiuXPn6rbbbpMk9ejRQ19//bWdJQAAAAAwlK3T3gICAhQdHS1JqqioUEpKyjV/63ZmZuZVlxcWFta7Pl9WWFiojIwMp8u4ZpX/b75YO3wPx5t5TB4bGBeA2nG8mafRws+OHTuUlJRUZVmXLl20ceNGlZSUKC4uTmVlZZo2bdo1bTc8PFyBgYFXLE9NTdX58+frVbMva926tSIiIpwu45qlpqZKkk/WDt/D8WYek8cGxgWgdhxv1yeXy1XtBZNGCz/Dhw/X8OHDr1heXFysGTNmKCQkROvWrVOzZs0aqwQAAAAAcHPkgQdhYWF69tlnFRAQYHfzAAAAAAxl6z0/hw4d0u7du9WtWzeNHTtWkhQaGqoXX3zRzjIAAAAAGMjW8NOzZ08dOXLEziYBAAAAQJID094AAAAAwAmEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABiB8AMAAADACIQfAAAAAEYg/AAAAAAwAuEHAAAAgBEIPwAAAACM4O90AQBQV+np6UpLS/N4/aysLElSQkKCR+tHRUUpMjLS4/YBAA2vPmNDfccFibHB1xB+ABgjODjY6RJgs/z8fOXl5WnXrl1Ol2KrvLw8NW3a1OkyAK/HuGCe6yb8mDrASQxyMEdkZCRn1wAAVTA24FpcN+EHAIDvCgkJUXl5ue655x6nS7HVrl27FBIS4nQZAOB1rpvwY+oAJzHIAQAAAHXB094AAAAAGIHwAwAAAMAIhB8AAAAARrhu7vmBM/jeFQDA5ZweFyTGBgDVI/zAUTxfHwBwOcYFAI2J8IN64dn6AIDLMS4A8Gbc8wMAAADACIQfAAAAAEYg/AAAAAAwAuEHAAAAgBEIPwAAAACMQPgBAAAAYATCDwAAAAAjEH4AAAAAGIHwAwAAAMAIhB8AAAAARnAs/Bw6dEjh4eFONQ8AAADAMI6En4sXL+qpp55SaWmpE80DAAAAMJAj4WfFihWaPHmyE00DAAAAMJTt4Wf37t26dOmShg0bZnfTAAAAAAzm31gb3rFjh5KSkqos69Kli4qKirRx40aPt5uZmXnV5YWFhR5v83pQWFiojIwMp8sAAK9i8tjAuAAAV2q08DN8+HANHz68yrLNmzfrv/7rv/Tggw+6l0VHR+vVV19Vq1at6rTd8PBwBQYGXrE8NTVV58+fr1/RPqx169aKiIhwugwA8Comjw2MCwBM5XK5qr1g0mjh52ruu+8+3Xfffe5/9+jRQ9u3b7ezBAAAAACG4nt+AAAAABjB0fBz5MgRJ5sHAAAAYBCu/AAAAAAwgq33/AAAYLe8vDzt2rXL9nYvXrwoSQoKCrK97by8PLVr1872dgHA211X4cfEAU5ikAOA6oSFhTnWdlZWliQ58vncrl07R/cdALzVdRN+TB3gKttlkAOAK8XGxjrWdkJCgiRp8eLFjtUAAKjqugk/DHAAAAAAasIDDwAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABiB8AMAAADACIQfAAAAAEYg/AAAAAAwAuEHAAAAgBEIPwAAAACMQPgBAAAAYATCDwAAAAAjEH4AAAAAGIHwAwAAAMAIhB8AAAAARiD8AAAAADAC4QcAAACAEQg/AAAAAIxA+AEAAABgBMIPAAAAACMQfgAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjGB7+MnJydEvf/lLjR07VjExMTp58qTdJQAAAAAwkO3hZ/78+RoyZIi2bdum6OhorV692u4SAAAAABjI387GcnNzdfjwYb388suSpHvvvVeDBg2yswQAAAAAhrI1/GRnZ+vmm29WYmKi9u3bp5tvvlmLFy++pm1kZmY2UnWeKywslCRlZGQ4XAkAwFswNgCA92m08LNjxw4lJSVVWRYWFqZDhw5pzpw5WrRokTZv3qy4uDi98sordd5ueHi4AgMDG7rceklNTZUkRUREOFwJAMBbMDYAgDNcLle1F0waLfwMHz5cw4cPr7LsH//4h8aNG6chQ4ZIkkaNGqVly5Y1VgkAAAAA4GbrAw86deqkDh06KC0tTZL08ccfq1evXnaWAAAAAMBQtt7zI0kpKSn6zW9+o6efflqtWrXSihUr7C4BAAAAgIFsDz9dunS5pnt8AAAAAKAh2P49PwAAAADgBMIPAAAAACPYPu3NW6Wnp7sfxHCtsrKyJEkJCQketx8VFaXIyEiP1wcANKz6jAsSYwMAeCPCTwMIDg52ugQAgJdhbAAA7+NnWZbldBF1UfllRd74JacAAAAAvENNuYF7fgAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABjB3+kC6sqyLElSSUmJw5UAAAAA8FaVeaEyP1zOZ8JPaWmpJOno0aMOVwIAAADA25WWlqp58+ZVlvlZV4tEXqiiokLFxcVq1qyZ/Pz8nC4HAAAAgBeyLEulpaVq2bKlmjSpepePz4QfAAAAAKgPHngAAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQAAAGAEwg8AAAAAIxB+AAAAABjB3+kCnPD+++/rhRdeUFlZmSzLUnR0tB555BFJ0qRJkxQaGqo1a9a43//cc89JkubMmaNJkybpzJkzatGihSSpqKhI//Ef/6HVq1frxhtvrNLO5e8tLy9XQECA5s6dq6ioKElSXFyc+vfvr/Hjx1d5r2VZsixLM2bM0IgRI+zokhrZ2V81bSsjI0OJiYkqKytTSEiIEhMT1bFjR8XFxSk7O1u///3v5efnJ0naunWrPvvsM61YsaLxOqYOnDjWvnv8FBcXa+HChfrqq68kSdOnT9fIkSMlSa+++qrefPNNWZYlPz8/TZkyRWPHjm3sbqmVtx9z33zzjZYuXaqjR49KkkJDQ7V48WLdcsstjdYntbGrzyTpnXfe0X//93+rvLxcTZo00bBhwzR9+nQ1bdrUhj1tHEVFRVqzZo3++te/qmnTprrhhhsUFxenXr166eTJkxo6dKheeuklDR482L3Oj3/8Y23atEmSNGzYMHXt2lWSVFFRoeLiYo0dO1a/+tWvqrRz8uTJKu+9dOmSevTooSVLlrj7ukePHjpy5EiV9/r5+am0tFShoaFKSkrSTTfdZEe31MjOPqtpW9///veVkpKiHTt2SJKioqI0f/58Sd/2ZUJCgu6//373epMmTdLs2bM1YMCAxumYWjhxrF3t+PGlzzfJN463ffv26ZlnntHFixdVXl6uqKgozZs3z7HPRrv6TJKKi4u1evVqffrppwoKClKrVq00Z84cDRo0yJ6dvRaWYc6cOWP96Ec/snJzcy3LsqyioiJr3Lhx1q5duyzLsqyHHnrI6tWrl/Xhhx+611m7dq21du1a9+t/+ctf3K+Vl5dbs2bNslatWnVFW99974EDB6z+/ftbx44dsyzLshYsWGBt2bLlqu89duyY1bdvX+vcuXMNtesesbu/atrWkCFDrL///e+WZVnW5s2brenTp1uW9W0/9urVy9q4caN7vS1btlgLFixokD7wlJPH2uXHzzPPPGOtWLHCsizLOnfunDV48GDr7Nmz1v79+63o6Gjr4sWL7teioqLcfewUXzjmlixZYq1fv969zrvvvmuNHTu2QfbfE3b22ZYtW6xRo0ZZWVlZlmVZVmFhoTVr1izrySefbLT9a2zl5eVWTEyMlZycbJWWllqWZVl79+61Bg0aZOXm5lrZ2dlWr169rCFDhliFhYXu9YYMGWJlZ2db2dnZ1pAhQ6ps88yZM1bfvn2tL774osry7763oqLCWr16tTVx4kT3su7du1/1vZZlWatXr7ZmzpzZMDteD3b3WU3b2rNnj/Xzn//ccrlcVklJiRUbG2vt3LnTsqxv+3LgwIHW6dOn3et993i3k5PHmmVVPX585fPNsnzjeHO5XNbgwYOtf/zjH5ZlWZbL5bKmTZtm/c///E9jdUuN7OyziooK66GHHrKWL19uuVwuy7Is6+DBg9bgwYMd+12riXHT3vLy8lRaWqpLly5Jklq2bKkVK1aoW7du7vfMmDFDS5cuVX5+fq3bu3DhgvLy8hQcHFzre3v37q3hw4dr8+bNtb63W7duatGihU6dOlXrexuT3f1V3bZKSko0d+5c3XbbbZK+PZv39ddfu19/+OGHtW7dOmVlZV3jHjYeJ4+1y4+f/v37a9KkSZKkdu3aKSQkROfOndPZs2dlWZYuXrzofm3t2rVq06aNB3vbcHzhmDt37pxcLpcqKiokSSNGjNCcOXOudVcbjJ19lpKSovj4eHXq1EmS1KpVKy1fvlx//OMfderUKW3dulVxcXHu90+aNEn79u3Tvn37NGHCBI0fP14LFizQ3r17NX78eI0fP15TpkxRbm5uPXvBc/v27VNOTo5+9atfyd//2wkRAwcOVFJSkvv/ODQ0VHfddZdWrlxZp21W/n61bNmyxvf5+flpzpw5OnbsmA4fPlzrdvv166cTJ07UqYbGZHef1bSt9u3bKy4uTgEBAWrWrJm6du2q06dPu1+PjY1VfHy8J7vZ4Jw81qR/Hz++9Pkm+cbxdvHiRRUVFbnH1ICAAC1atEj9+/f3dLfrxc4+++yzz3T69Gk9+eSTCggIkCT17NlTM2bM0PPPPy/p32OB9O3VtR//+MeSvp0FNX36dA0fPlwfffSRVq5cqTFjxmjcuHFKSUmpf0dchXHT3m677TYNHTpU99xzj26//XYNGDBAo0ePVlhYmPs9/fr1U35+vhISEqpME6kUHx+voKAg5ebmKjg4WCNGjNDkyZPr1P6tt96qTz75pNb3/elPf1J5ebn7cqNT7O6v6rYVEBCg6OhoSd9eek1JSdE999zjfj0sLEzTp0/XwoUL9fvf/76B9r5+nDzWLj9+Lv+QSk1NVUlJibp166YuXbpo69atuvvuu3XHHXdowIABio6OVocOHRpk/z3lC8fcjBkzNGvWLP3hD3/QwIEDNXjwYI0ZM6YBe+Ha2NVnubm5OnXqlPr06VNleXBwsLp166aDBw/WWOeJEyf08ccfq3Xr1po0aZJ++9vfqk+fPtq0aZMOHTqkH/7wh553Qj0cOnRIvXv3VpMmVc8HVk5RPnnypKRvB+nRo0drz549VaaJSFJOTo6io6PlcrmUl5en3r17KyUlpU7T0wICAhQWFqavvvrK/cfo1ZSWlmrHjh268847r3UXG5wTfVbdtm699Vb3zydOnFBqaqpef/1197KpU6fqww8/1Jtvvlll+psTnDzWLj9+fOnzTfKN4y04OFjTpk3T+PHj1blzZw0YMEDDhg1Tv379GqQPrpWdffb5558rPDzcfetBpR/84AdXHW++KyQkROvXr9epU6e0Zs0avffee3K5XFq0aJFcLpcCAwOvef9rYtyVH0launSpPvroI02cOFGnT5/W/fffr507d1Z5z2OPPabPP/9cu3btumL9ZcuW6Z133tHatWtVUFCgqKgod9KtjZ+fn5o3b37V1+Lj4xUdHa1Ro0Zp/fr1+t3vflenMzmNze7+qmlbJSUlevzxx1VWVqZp06ZVeS02NlaWZbnnqnoDO/uutuNnx44dSkxM1Nq1a+Xv76+AgAA9//zzeu+99zR8+HAdPHhQY8aM0f79+xu0Dzzh7cdceHi4du/erbVr1+qWW27RSy+9pAceeEBlZWX13HPP2dln5eXlVywrLS2ttcbOnTurdevWkqShQ4dq9uzZeuqpp9S1a1fHgo8kNWnSRJZl1fq+Vq1aKSEhQYsXL1ZRUVGV10JDQ7V9+3alpqYqOjpapaWlGjhwYJ1rqG5sqPzjIzo6WmPGjJFlWZo3b16dt9tYnOizmrYlSceOHdMvfvELLViwoMr9Kf7+/lqxYoWSk5OrzBhwgt39Vtvx4yufb75yvM2YMUPp6emaNm2aiouLNXXqVG3cuPGa9rWh2Nlnfn5+Ho8Lktwn1Dp06KDAwEDFxMRo48aNevTRRxs8+EgGhp9PPvlEqamp6tChg+69914lJycrPj5eb731VpX3BQUFKTExUUuXLlVBQcFVt3XnnXdq0qRJWrBgQZ0/FI4cOVLt1Zxly5Zp+/bt+uMf/6hXX33VK24Sc6K/qttWcXGxHnnkEZWVlWndunVq1qxZlfWaNGmipKQkr5n+Znff1XT8vPLKK1q5cqU2bNjgPrO8bds27d27V2FhYXrwwQe1fv16/ed//qe2b9/eQD3gGW8/5izL0m9+8xuVl5erf//+evTRR/XOO+8oLy9Phw4daphOuEZ29Vnbtm3VqVMn/e1vf6uyPDc3V9nZ2e4zf5cPuJcPfpf/cT958mS98sor6tSpk55++mmtW7fO4/2vr/DwcB06dOiKPxSeeeYZ/eUvf6my7Ic//GGN00SaNGmi+fPn6/z583rppZfq1H5JSYmOHz9eZZpipco/PrZv364dO3Zo5cqVCgkJqduONSKn+qy6bWVkZGjy5MmaN2+exo0bd8V63bt394rpb3b3W03Hj698vkm+cbzt379fr776qtq2batRo0YpKSlJKSkpdbrVoTHY2Wd9+/ZVZmbmFWFn//796t27t/vflbV8d2ypHBv8/f21efNmzZ07V/n5+YqJidHx48fruMd1Z1z4ad68udasWeO+3GdZlr744gvdfvvtV7y3X79+GjZsWJXL5981ZcoUXbx4scb3VDpw4IA++OADTZgwwfMdsJlT/XW1bT3xxBMKCwvTs88+W+0Z6crpbxs2bKjL7jUqJ4+1y+3atUsbN27Ua6+9ph49eriXl5eXa82aNe57LcrKynT8+HH17Nnzmrbf0Lz9mPPz89OXX36pDRs2uOdN5+TkqKyszH0fjN3s7LNHH31UiYmJys7OlvTtH1Dx8fEaMWKEbr75ZrVp00ZffvmlLMtSdna2jhw5ctU27rvvPhUXF2vy5MmaPHmyo39Y9evXT+3atVNKSor77OWf/vQnbd269aqBJC4uTp9++qlycnKuuj1/f3/Nnz9f69ev19mzZ2tsu6KiQs8995z69u3r2PHjCSf77Lvb+vrrrzVr1iytXr3a/STLq5k6dary8vKuCO92crLfvstXPt8k3zjegoODlZKSUuXevWPHjl31c9gOdvZZv3791K1bNyUmJroDUGZmptatW6eZM2dKktq0aaMvvvhCkq46+0D6dqreQw89pB/84AdasGCBunbt2ijhx7h7fgYOHKjZs2dr+vTp7v+gu+++W7Nmzbrq+x977DGlpaVVu72AgAD3HwPR0dHuKR2V4uPj1aJFC/n5+SkoKEjJycn6/ve/33A71Mjs7q/qtnXo0CHt3r1b3bp1cz+KOTQ0VC+++OIV68XGxl4x3ccJTvbd5dauXSuXy6Xp06e7ly1btkz33nuv8vLyNHHiRPec4JEjRzoezn3hmHvmmWeUlJSkoUOHKigoSK1bt9aaNWscOyNvZ5+NHDlSTZs21dy5c1VSUqLy8nKNHDnSfXzddddd2rJli4YNG6bOnTsrIiKi2hri4uLk7++vwMBALV261NPdrzc/Pz89//zzSkpK0qhRo+Tv7682bdrohRde0I033ugOlZUqp4k8/PDD1W4zMjJSd9xxh373u99p+fLlVV6rnIokfRt+br/99jrNi/cmdvdZTdvasGGDXC5Xla82iImJ0cSJE6usVzn9bfz48Z7scoNwst8u50ufb5LvHG8rVqzQwoULVVRUJD8/P/Xt21dLliyp5957xu4+S0lJUXJyskaNGqWmTZsqODhYTz/9tPuR8o888oji4uK0ZcsWDR069Krb79mzp+644w6NGjVKQUFBuv322xUZGVnPnriSn1WXCYEAAAAA4OOMm/YGAAAAwEyEHwAAAABGIPwAAAAAMALhBwAAAIARCD8AAAAAjED4AQD4jGeffVbbtm2r8T1bt251f1v9d02aNEnvv/9+I1QGAPAFxn3PDwDAd82dO9fpEgAAPozwAwCwzbx589SzZ0/3F+m99tpr2rt3r0JDQ/V///d/Ki4ulmVZWrZsmSIiIhQXF6f8/HxlZ2frRz/6kc6fP69bb71VDz/8sN566y298cYbKi0tVUFBgaZOnaoHHnhAknT27Fk9/PDDysnJUceOHZWQkKD27dtXqeV///d/tXr1al28eFF+fn6aM2eOhgwZYnufAADsw7Q3AIBt7rvvvirT1rZu3aoePXooJydHb7zxhlJTUzVu3Di9+OKL7vdcunRJ7733np544gn3suLiYm3evFkvvPCCtm3bpuTkZD399NPu148fP64lS5bo3XffVffu3a/4NvKCggI9+eSTWrVqld5++22tW7dOv/3tb3X69OnG23kAgOO48gMAsM2AAQPkcrn0+eefKygoSLm5uZo5c6aOHz+u119/XdnZ2dq3b59atmzpXiciIuKK7bRs2VLr169XWlqaTpw4ocOHD+vChQvu1++66y6FhYVJkiZMmKAJEyZUWX///v06e/asZs2a5V7m5+enI0eO6Oabb27o3QYAeAnCDwDANn5+fpowYYK2b9+uZs2aacKECUpLS9Py5cs1ZcoUDR06VF26dNE777zjXqdFixZXbOfMmTP6+c9/rvvvv18REREaNmyYPv74Y/frTZs2df9sWZb8/asOd+Xl5eratas2b97sXvbPf/5Tbdu2bcjdBQB4Gaa9AQBsNW7cOH300Uf64IMPNH78eO3Zs0dDhgzRAw88oN69e2vXrl0qLy+vcRuZmZlq27atZs6cqbvvvtsdfCrX27dvn3sK22uvvabIyMgq699xxx3KysrSX//6V0nS3//+d/3sZz9TTk5OQ+8uAMCLcOUHAGCr9u3bq2fPniorK1OHDh0UExOjxx9/XKNHj1bTpk3Vr18/7dy5UxUVFdVuY/DgwXrrrbc0bNgwBQUFqU+fPmrbtq2ysrIkSd27d9fChQt17tw5denSRU899VSV9du2bau1a9dq1apVcrlcsixLq1atUseOHRt13wEAzvKzLMtyuggAAAAAaGxMewMAAABgBMIPAAAAACMQfgAAAAAYgfADAAAAwAiEHwAAAABGIPwAAAAAMALhBwAAAIAR/h99zfaJErBRIgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "\n",
    "df_diff = pd.DataFrame({'SNR DIP':SNRs['dip'],\n",
    "                        'SNR N2N':SNRs['n2n'],\n",
    "                        'SNR P2S':SNRs['p2s'],\n",
    "                        'SNR N2S':SNRs['n2s'],\n",
    "                        'SNR Ours': SNRs['our'],\n",
    "                        'CNR DIP':CNRs['dip'],\n",
    "                        'CNR N2N':CNRs['n2n'],\n",
    "                        'CNR P2S':CNRs['p2s'],\n",
    "                        'CNR N2S':CNRs['n2s'],\n",
    "                        'CNR Ours': CNRs['our']})\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set(rc={'figure.figsize':(14,5)})\n",
    "sns.set(style=\"whitegrid\")\n",
    "set2 = sns.color_palette('Set2')\n",
    "my_pal = {\"SNR DIP\": set2[-1], \"SNR N2N\": set2[0], \"SNR P2S\": set2[2], \"SNR N2S\": set2[3], \"SNR Ours\": set2[1],\n",
    "          \"CNR DIP\": set2[-1], \"CNR N2N\": set2[0], \"CNR P2S\": set2[2], \"CNR N2S\": set2[3], \"CNR Ours\": set2[1]}\n",
    "\n",
    "ax = sns.boxplot( x=\"variable\", y=\"value\", data=pd.melt(df_diff), fliersize=0, sym='', palette=my_pal)\n",
    "\n",
    "# plt.savefig('quantitative.svg', dpi=500, bbox_inches='tight')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "754147cbc932618733e4b5b22f4733dc08362733592b3e8aa661889d22935194"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 64-bit ('tsne': conda)",
   "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.6.10"
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}