{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ivim\n",
    "import os\n",
    "import tempfile\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Define output folder (can be any valid folder)\n",
    "folder = os.path.join(tempfile.gettempdir(),'ivim_example')\n",
    "if not os.path.exists(folder):\n",
    "    os.mkdir(folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Protocol optimization\n",
    "Before acquiring data we might want to optimize our acquisition scheme. This can be done using a CRLB-based optimization. The input needed is expected IVIM parameter values of the particular IVIM model (referred to in terms of a specific regime, here the usual biexponential model with D* corresponding to the diffusive regime). We might also want to set a maximum b-value to avoid low SNR or kurtosis effects at higher b-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b-values:    [  0.  44. 222. 800.] s/mm2\n",
      "Proportions: [16. 30. 36. 18.] %\n"
     ]
    }
   ],
   "source": [
    "D = [0.8e-3, 1.0e-3]\n",
    "f = [0.1, 0.15]\n",
    "Dstar = [20e-3, 15e-3]\n",
    "regime = ivim.models.DIFFUSIVE_REGIME\n",
    "bmax = 800\n",
    "res = ivim.optimize.crlb(D, f, regime, bmax=bmax, Dstar=Dstar)\n",
    "bopt = np.round(res[0])\n",
    "aopt = res[1]\n",
    "print(f'b-values:    {bopt} s/mm2')\n",
    "print(f'Proportions: {np.round(aopt*100)} %')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is a tuple with two vectors corresponding to the optimized b-values and the fractional time one should spend acquiring data at each b-value. In other words, let's say that we have time to do 15 acquisitions, then we should use the following acquisition scheme:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acquisitions at each b-value\n",
      "  0: 2\n",
      " 44: 5\n",
      "222: 5\n",
      "800: 3\n",
      "\n",
      "Resulting b-value scheme: [  0.   0.  44.  44.  44.  44.  44. 222. 222. 222. 222. 222. 800. 800.\n",
      " 800.] s/mm2\n"
     ]
    }
   ],
   "source": [
    "n = 15\n",
    "print('Acquisitions at each b-value')\n",
    "for ai,bi in zip(aopt,bopt):\n",
    "    print(f'{int(bi):3d}: {int(np.round(ai*n))}')\n",
    "\n",
    "b = np.repeat(bopt.T,np.round(aopt*n).astype(int))\n",
    "print(f'\\nResulting b-value scheme: {b} s/mm2')\n",
    "\n",
    "# Write the result to a bval file\n",
    "bval_file = os.path.join(folder,'example.bval')\n",
    "ivim.io.base.write_bval(bval_file, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "We can now go ahead and acquire our data. Here we will do so by simulations. Based on available maps of the IVIM parameters we can create noisy images based on the specific IVIM model. A simple square digital reference object is used here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sz = [20,20,100]\n",
    "pars = ['D','f','Dstar']\n",
    "files_sim = {}\n",
    "for par, vals in zip(pars, [D, f, Dstar]):\n",
    "    im = vals[0] * np.ones(sz)\n",
    "    im[sz[0]//4:3*sz[0]//4,sz[1]//4:3*sz[1]//4,:] = vals[1]\n",
    "    files_sim[par] = os.path.join(folder, f'sim_{par}.nii.gz')\n",
    "    ivim.io.base.write_im(files_sim[par], im)\n",
    "\n",
    "noise_sigma = 1/40 # SNR = 40, given S0 = 1 as default\n",
    "outbase_sim = os.path.join(folder,'ivim_sim')\n",
    "ivim.sim.noise(files_sim['D'], files_sim['f'], regime, bval_file, noise_sigma, outbase_sim, Dstar_file = files_sim['Dstar'])\n",
    "im_file = outbase_sim + '.nii.gz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAGBCAYAAAAZl6lgAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAUBpJREFUeJztnX2MV9Wd/9+MPIsMTwNjKWQ1Bh+yq3V9wJSOQ9WC2oJui4ghKnHbNa4mhga3pqtkE9MY21g1aWKoDxhlZcOyoNtJ3dk1OkiFrYmFtcQtshSr2zDOwDAwIk/O3P3D38xvvu/7nu/nfB8Gxvp+JSbc7/d77z3nfM45c7zv9/2cYVmWZTDGGGPMF5qaU10AY4wxxpx6vCAwxhhjjBcExhhjjPGCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDEAhp/qAvwp8Nprr+Gtt97Ce++9h1WrVmHEiBGnukgGjsvnBcdp6OLYDF0GIzaD8oSgpaUFK1asyH3e1NSE559/Hl1dXZg9ezbGjRuHHTt2AAA+/PBDLF++vKL7rlixAg0NDVi6dCmOHz8e/r6rqwtLly6t6J4AcNVVV+H+++9HbW0tjh49WvH1BosoLm+//TYaGhrQ2NiIxYsX48SJE47LIBLFY8eOHZgzZw4aGxvxzW9+Ex9//HFV4gEAa9euRV1dXfjZQFQrRsDQjFMUm176t5ljc3KIYvP++++jrq4Oc+fOxdy5c9He3u7YJHJSJYOnn34aS5YswZgxY9DU1IRFixb1fTdjxgy0trais7OzrGtv27YNra2t2Lx5My644AKsX78+POfVV1/F1VdfXdb9mKeeegrz58/HGWecUZXrnUx64zJ9+nQ0Nzdj06ZNOOecc/DSSy85LqeA3nice+65ePPNN7Fp0yZcfvnl2LhxY8XxAICenh6sX78eM2bMKPpZMaoZI+DzE6fe2AD5NnNsTi39Y9PY2IiWlha0tLSgrq7OsUlk0BYEO3fuxMKFC3HZZZfhnXfeQWdnJ44cOYKRI0di+PDhckXV0NCA5ubmsu63detWzJs3DwBw7bXXYsuWLbnvZ8+ejcbGRqxcuRIA0NzcjPnz56OlpQXz5s3DDTfcgIsuugjr16/vK3tbW1v4/bPPPoumpibs2LEDHR0dZZX/ZFEsLvX19Rg7diwAYMSIERg+/DNFyXEZPIrFo/8jwE8++QTnnXcegMriAQAvvvgiFi1ahJqamqKf9TKYMQIwZONULDaAbjPH5uQQxebNN99EQ0MDfvjDH6J3Q1/HJoFsEHj99dezOXPmZD09PdnOnTuzhQsXZr/+9a+z7373uwW/u/3227Pf/va3fcevvPJKtnLlyoLftLW1ZY2Njbn/9u/fX/C7H/3oR9nGjRuzLMuyXbt2ZbfcckvB9w888ED2i1/8IsuyLOvu7s6yLMsWLlzYV95rrrkmy7Ise/rpp7Mbb7wxy7Ise/zxx7NVq1aF339eSI3LH/7wh+yrX/1qdvz48SzLHJfBIiUe//7v/5595StfyS677LK+tq0kHp9++mm2YMGCrLu7O7vkkksG/Kw/X8QYRbEZqM0cm8Enis3Ro0ezjz/+OOvp6cn++q//OvuXf/mXLMscmxQGzVR48cUXY9iwYZg1a1bfimb06NHR4iT3WV1dHVpaWsL7TZw4EYcOHQIAdHZ2YtKkSQXf33333Xj44Yexbt06LFmyBGeddVbf/3EBwIUXXggAmD59esG/d+/enfT954UoLocOHcKtt96K1atX9/0fquMyeETx+MY3voFt27bhxz/+MX7+85/j/vvvrygea9asweLFiwv+j0Z91p8vaoyKxWagNnNsTg7FYjNq1CiMGjUKAPCd73wHW7duxbe//W3HJoFBWxBs374dWZZh9+7dmDp1KmbNmoU9e/YUPWfPnj04//zzCz5rb2/HTTfdlPvthg0bCv64XHHFFXj00Udx2223obm5GXPmzCn4fW1tLZ544gkcP34cl1xyCe644w5cd911fd8PGzZM/ru3E0Xff14oFpfu7m4sXboUK1euxKxZs/rOcVwGj2LxOHbsWN/EVltb22fIrCQe7777LrZt24Y1a9Zg165dWL58OUaOHJn77LHHHus754sao2KxUe342GOPOTYniWKx6erq6tPU33jjjb54ODYxg7YgqK2txYIFC/DRRx/hmWeewYQJE1BTU4OjR49i9OjRuP7667F9+3bs3LkTd955J5YtW4Y33ngDTz75ZMF1UldwF198Merr69HQ0ICZM2fivvvuK/h+1apV2LBhAw4fPoxly5Zhy5YtuOeee6pZ5c8FxeKyceNGbNmyBV1dXXjooYdw11134eabb3ZcBpFi8Xj11Vfxk5/8BDU1Nairq8Nzzz0HABXF45FHHun796WXXlowgQ302Rc1RsViM1A7OjYnh2Kx+dWvfoUHHngAY8eOxVlnnYWHHnoIgGOTxKCJEYKmpqZs9erV8rsPPvggu/fee09aWdauXXvS7jXUcVyGFkMpHoovcowcm6GLY1M5w7JsiD8bMsYYY8yg49TFxhhjjPGCwBhjjDFeEBhjjDEGXhAYY4wxBl4QGGOMMQZeEBhjjDEGXhAYY4wxBl4QGGOMMQYlpC7+/e9/X3Dc09OT+013d3fB8eHDhwuOx40bV/R7APj0008LjntzuQ90DQB9m+f0wntDc7kA5DakOHHiRNFjBd93woQJBcf9t6/tpTcffS9/8Rd/Ed4n4r/+678KjseMGZP7DZelq6ur4DhlkxbenIr3Fh8/fnxY1hR6Nyvppb6+vuD46NGjBceqrNxPOFbcz4B8vzn33HPjwhbhV7/6Vfgb1Tf7w3VTdeUxcuTIkYLj3i1h+3PaaacVPeY2Tikb30fNEVzWTz75JLwPt9FVV10VnhPBsVFl5THBfYbHGY9tIN+uPO7UPvaqb/anf177Xo4dO1ZwzLFI6Ud8XS479ysAfVuk9/L1r39dlDgd3p5YbYjHceFyczlVG3PduM3379+fO4fnd+7LKm58H+4jqt8x3P95Lud7APlxdfXVV4f3AfyEwBhjjDHwgsAYY4wxKEEy+PjjjwuO1WMrfjTDj935MZXaR3rKlClF78uPxtR1+HGZOocfK/FjGX7krB6h8iMkfoTIZQc+26Wr2nAsVFn5cRY/duJj9SiX2/n0008vOObHh0BeFuI+oB531dXVFRxzbPi+KY/qosehgO6PlcBtOHbs2Nxv+JEo9xmuG9ddXYP7u5Lm+HEnP7rk71VZuA15vKtH6HwfjovqQ9WOC5CPjbovtzX3Xb6GerzNfZfHvxpnPH/xsYonl5/bjMe3miP4N9wX1eNtJU9WQsojdO6H3O+4LdQ8HMnKSmbg3/A4U39n1Jjvj5LzGK5fShup8ZuCnxAYY4wxxgsCY4wxxnhBYIwxxhh4QWCMMcYYlGAqTDH2sEEoMrIpwwWfwwYRZWZkYwa/x67yAXB92IQRGdmAvFmFjYqqfsrgUils7FHmlugdajaqqDwMfF1uV/WeMhut2JilTElsRuNzOjo6ipZDXYP7jerPqt0qgY1JyvzI7RwZN9X743wNbg9lZGUzW2SYAvL9ma/BZVU5Fng8p7yXze/uVwPud2peYfMezzNcF1VfnlcOHjxYcKyMl9E11H04Nlwf7hMpfT0yxKnrVgrPFyoufE+OA3+vxh33Mz5H9UO+D8dy4sSJuXP43lFuC9WeHIfIeKzKmoqfEBhjjDHGCwJjjDHGeEFgjDHGGJTgIYgSqAD5ZBesubHGrJLDRNpeiobKmptK/qD07v5w/ZS2wxox63IpCSSqAZdV6eNRrnpuI7VnBOtS0b4TQL7dOFbt7e25c6KkGlFSJSDO/608EuUm8xgIvodKBsPtzG3KfUjFNkruo+7L2jUfq/hzP+M25jGl9F9uEx7vvOcEoOeJSuGypSQIGgxfSopmzO2s7sNtH/kbUvZ34XZXdUnxQJQC+11SvEz8NyHyIKlrcCzV3M2/4cR0iug+3MfUWOVzuE3UWC03Ln5CYIwxxhgvCIwxxhjjBYExxhhjUIKHgDUJpWPxZ6x1sF6itEF+95d/ozQ3/oz1MaUPqo1I+sPvlCrNjXUmLqtqo8HwFZTzzmmkVav2iTZASokn9wGV3yEqK7er8oNEm7soqp2HgNtQ6XrRPvYp7ylHOQRScndwH1JtwVrspEmTCo5TdEuuT5TbAojHajnwfdU79zx/jR8/vuA48mEoIg0ZyI8r9iqp9uDxG/kdUnKzcDlUzhDVtyqB54uUPATcL3k+UOXm9lKbGTGs73NfVvkAOL4cu87OzoJj9XeG58gUX536LAU/ITDGGGOMFwTGGGOM8YLAGGOMMSjBQxC9tw7kc46z5sbvebNeBOT1Hj5W57BewsdKy2fNlK/LGpPSVPmdWdZ2VK4GpWdVCrer8ilEubpZL1UaFNcn5d1dbmfuN+qc6J1p1tRU/oBy3hlP2a+jFPieymMReRv4HNVeXFdu4wMHDuTOYS0zJd8Bv+8ceRWU1s3zCB8rLb/aOnXqNXkMKH23Pyq+Ue4VNa/w++2cm0F5Fbhs7O/g8Z3iieJ48x4xQPVzdzApvgweQ9z/U3IZ8PyvvE18zrRp0wqOVftEORP4+5RcFil7yER9dSD8hMAYY4wxXhAYY4wxxgsCY4wxxsALAmOMMcagBFNhSkKZaEON/fv3F70mkDfd7Nu3r+g1gbwRia+h7hOZmaKkJADQ0dFRcMwmKmUQUZu3VAonEFFl5fqyWYcTcygjDtcn2gwFiI2Hqo2ipBopiYm4L/J9lZlP1bkSuB7qnvybKBGRapsoAdKUKVNy53Bd2TCqDIHcd7lsfA6PD/UbjpMy+w3GmIn6MpBva64Pjxm1MQ2biLmd1TzKJlBuk5RNxKINfpS5mY11Bw8eLDieOnVq7hxV50pgA51qH74nj3+e79XfjGgTPfU3g9u4tbW14FjNf1HiNTaQqjkoSiKn5oSURGwKPyEwxhhjjBcExhhjjPGCwBhjjDEowUPA+ggnIQLihBDlJLFgLUTpNPxZShIO/oy1TNZtlE7NGyDxOWpTpZRNNEqFNSR1X64feyaiTYiAfLtyu6dsRMJlU9ot63tcFtaUVTIbLhtrt0p3K2eTqGKkbJYUbTrD7cOboQD5+rPGqJJ5sbbJ12BtU92b48L3UXqoilV/lCadsgFWqaQk1eI+wvMX9xfVf9g3xX07ZT5jVMIznou4HbkfKQ9B9BsVm2onJkrR1KO/K1xONf/zXM2xU/MS3zcl2RnDfzdT/h5ECc6U36fcZF5+QmCMMcYYLwiMMcYY4wWBMcYYY1CCh4D1E6VbsC7D2l/0PiWQ13/4Pko/Y32f76vOYe2Gr8FarnqvkzVVvqbS6bgNqkHKO6f8XjLrbpF2DeTLHmmq6r7Km8BwLMrZIInLyrFRZa12bFjrU/eMfAb8vcoxoepfyvdAvk2VD4W1WG4v9gfw5jhA3mcQeSiAtPKXCl9Taa5qjutPNN8B+Zhzu6pcDRxj9jKofhTNm1w/5QeINoVT46PaseH5Qt2Ty8V9l+uq8lhEuQpUvaLN+VTf5fpw2flvk7pGNL+rNopyuQyEnxAYY4wxxgsCY4wxxnhBYIwxxhh4QWCMMcYYlGAqjDbPAPImOjZdsMlKJYxgYw6bcFQiBzZmpGy4weXnc9iYoowbXP6UzWxSkleUCpvwamtrw9+wSZSNSiq+nDSEY6MMctxG5Rjt2AAUbVyk4IQgKjblJvMYCO6XylDJRrwooUxKQiWOnUp0w/fh/qAMgRxLHt/c79ra2nLXiOKvEhcNhhGXy5EyLrmNeC5SBkGOBRsGVT/kccX9Uo2ZKKkWj39leOb6cLurpFjTpk0ret9S4XuquETmVm6vlMQ9KXMKb/Y0adKkgmM1/0Xm9CjZmzqH+4caM+WaPf2EwBhjjDFeEBhjjDHGCwJjjDHGoAQPQYqmHiVZiLQfIK+PcDIItdkF6z+slyltmMvCWtWBAwcKjtVmL6zd8DXUhkjR5i7lECXmUL+JdFmlVXObpGw8FelfinI2wGEibU5pbNVOssL1UH03ShgTJWEB8rGKEgYBeW2b9f+UzX74ujym6urqctfg+PM1VWyr7e2o1jVZU1btHCXRUu2csgESw/4Fvm9KX4zqo+qXMp5LIaV9eE7h33DfTkkyxHOb8mFNnjy54Jg3rlJ/I/g3TEqsOVbss1NePLX5YAp+QmCMMcYYLwiMMcYY4wWBMcYYY1CChyB61xPI6y6sSaVsbMP3YU1Fad/Ru83qHV3Ww6INktQ7pnxd1kOVtqPe5a2UlHf7eYMPrh9r1UrLUrkJ+qM8Exxz1v/UZj3qfe7+cP1YD1S/ifqVKlulpGzsxBosvzPNfUjVlfVCvobqc6x3cvuoPhT5UFI20GE9lM9R3pVqxwXIxz9Fp4+0alV2rm+0qQ6QzwfA+rbS/3k8Rzq6im8016o4RPkPKkVdP8rTwLq8mpd4Pmd/gLovz0spfTcaMxw3FVtu9yiHCKDzTKTgJwTGGGOM8YLAGGOMMV4QGGOMMQZeEBhjjDEGJZgK2fylTHZs3og2WVFGnsiYp8wSbLJIMQSyMSUyPKZsGJSSMCfFvFQqbDpRSSl4sxo2KqUkyODfsPlJJSlhcw7/RiXu4LKyaYYNkCq+UTnUOSmm11JgA5HamItjxWOG+3JKopYUeJzxMZtQgXx92ODIx6pcUXIydY4yTVVKymY/UXIrNm+mGFXZiKYMb9wHuGyqrLzRFG+8w+MuZYMkVTYm5TelkGLU4/pzv+O6qrhwP+T5QPX/KHmXMryzeZ1/wxuA1dfX567BsY02RAPK30TPTwiMMcYY4wWBMcYYY7wgMMYYYwxK8BCwrqX0EtZyIg+B0oeiJDspm9KkaLd79+4tOObEFCn6EJc/JUnHYCTyYP1LJY1iTYl1ei6X0ilZ7+bNrJRWx+0YbZCjiDbWUsld2A+Q4g8oV3cbCK4ba4FAXoeOEpmkJC6JNu4C8mNTlY2J9FA+Tqkvl1WNj8HY3CjFhxL1B9ZuVdl5nLHmrvo/twkfK58Fz73c9inepSjRlJp7qz2flVNOriuXW8WW2zCaH9V1UsYZ14d/w8fKh8VzZpQgTt03FT8hMMYYY4wXBMYYY4zxgsAYY4wxKMFDwO+Lp2zUEr0fr7T9aJMVtVEL50jga6j3NNW9i923vb0995u6urqCY9Z3lb6mfBOVwpqZ0lxZ/4ramfMUAHlfAWtoKv8B94GUzV24HaONiZRexvVlrS6ljSolxcfA+mf0TrXaDCrSGJW2ydfl99bV++XczpHGrDZi4s+ifAjA4Phu+D4pm3lx/+ccEaov8/vsHAtVX/5N5MVSRBtN8dwFxDkFqrmJzkCwL0mN0ylTphQc89+mqJ8C+f7Pfw/K2aiNyw7kPSI8x7BXgT02qqwcS9Xvyv074ycExhhjjPGCwBhjjDFeEBhjjDEGJXgIVC5/JtKyWddIyWXAGpx6b5e1LdZ6lPbFOgxrN6xbqvuytseaktJQVY7sSlEaEsP14WPWLVV9WS9kzWzixIm5c6I9I1L2f4i8C6ofsTaX4oFR8aoE9gekvD/O/YPbVLVXpJmquLCWqbRLhnVVHleswyp9ubW1teCY46R8HKzVDwbqHfsoDwGfo7Rqjg33gRTPBMemnHwZPFaVRh75itR9q73PBNc1RVOP9hg4ePBg7ho8JnguZ58CEM/dapxFfo+UuYyvweNM9d2UvCIKPyEwxhhjjBcExhhjjPGCwBhjjDHwgsAYY4wxKMFUyIYZZTBhAwgbWdhAocwQTIqBhk0XbIhTBho287BhhMuaYmziZD6c7AVIq3OpsMlEbZARmZm4XCnJfvg3ymAUbUSSYuTj+7IhUMWXy8KxUOa1FGNdKXDdlPmTy84mO66HGg/RxiyqfbjPcFIp1U/ZJMhtyOeo5Fa8iRibuSZMmJA7RxmtKoUNcymJ1vi4nI25uL4q6RLXl5NRKWMptz33NY6d2iCn1PgC1Td8cnuo2Ecb7XE5U0y1jGrjqB+m/E3kcZay4R/3zZQEauVu1OYnBMYYY4zxgsAYY4wxXhAYY4wxBiV4CFI27uGNh1g/S9nYhjVm1rVUshM+J2WDDKUZ9of1QaW5sU7DZVNaforuWCqceEO1EcePyxElEALy7ZySAIo1MtbUVB+IfBZcVuVD4OuyRyRFM66UaLMYIK/lR4maVF/mfljOBlJ8XdVPOS78mygJD5DvQ3wNtXnZYCQm4jGi9OGoTbjfcfIrIG4T1Xf5N7yxmvL38Dk8X3G/UnNENK6Udq3mxUrgcqqEeBwr7v/RpnpAPpY8RlKSTKWcw3+LovirOYOvwXOGmsvKTYDnJwTGGGOM8YLAGGOMMV4QGGOMMQYleAgYtXkCv0PJ2h9/r965ZG2Lr6H0oOg9XaVzsZbD9+WyKW23rq6u4Lijo6PgWGlK1dbcgHxdlGbMn3E5+HtVX74P+0iUh4DbgPuNekeYdWTWrll3U+3MeidrkazLAfn39ysl2vwKyPeZaHMXpeOytsmeEpVfgePLbaz6abTZDWu7Ki5clii2wOBsCMbtqMaMild/Il8OkJ9XUjad4fhxm6h3zFl75rbneVT1I74PjxEVm2rnVeG5XJWT+6H6O9IfpbFzvDnWql78G76u6u98Dsefxy7PB0B+7uJrqD5Vbk4VPyEwxhhjjBcExhhjjPGCwBhjjDHwgsAYY4wxKMFUGBmK1GdsCGGDhUoGwmYONpix6QTIm3DYQJNi5GODUGQ6UWVjA5EypijzUqWw2Ucl2OFkFpzsZNq0aQXHbW1tuWtwu3L8lKmQ+w0f79u3L3cOm7O4XdnspBK1RImluD2AtE1DSiFK3APkxwybirgPKYNUtLmXMuKqhC/9UUY+Nm+xCZP7nUqywn2V+4MyZg1GMi8eh6rv8m8io5aaZ3iMcF3UHMjtyvFTJjruz9GGP6ovstGUyzEYm0wx3OaRYRCIE6Spa0Tzu5qneaxyX1Vtyu3O1+DvlbE5Sm6VkqwsFT8hMMYYY4wXBMYYY4zxgsAYY4wxAIZlKSKNMcYYY/6k8RMCY4wxxnhBYIwxxhgvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDLwiMMcYYAy8IjDHGGAMvCIwxxhgDYPipLsCfGq+99hreeustvPfee1i1ahVGjBhxqotk/h+OzdDHMRqaOC5Dn2rEqGpPCFpaWrBixYrc501NTXj++ef7fnP11VejsbERL7/8Mj788EMsX7684nuvXbsWdXV14WcD0dXVhaVLl1ZcDgC46qqrcP/996O2thZHjx6tyjUrJYrNW2+9hblz52Lu3Lk499xzsXz5csdmEIni0dPTg9tvvx0NDQ248sorsXv37orj8fbbb6OhoQGNjY1YvHgxTpw4kfxZCtWK06mMURSX7u5uLF26FHPnzsWyZctw4sQJx2UQieLR1dWF2bNnY9y4cdixY0ff9ytWrEBDQwOWLl2K48ePA0BFcSolHuXEaSjFaNAlg6effhpLlizB0aNH8eijj+KVV17Bpk2bcMMNN2DGjBlobW1FZ2dn2dfv6enB+vXrMWPGjKKfFePVV1/F1VdfXXYZmKeeegrz58/HGWecUbVrDga9sbn88svR0tKClpYWNDQ04MYbb3RsTgG98di+fTuOHTuGzZs348EHH8TPfvaziuMxffp0NDc3Y9OmTTjnnHPw0ksvJX+WQjXjNNRi1BuXDRs24Oyzz0ZLSwsuuOACbNiwwXE5BfTGY8yYMWhqasKiRYv6vtu2bRtaW1uxefNmXHDBBVi/fj0AVBSnUuJRTpyGUoyquiDYuXMnFi5ciMsuuwzvvPMOOjs7ceTIEYwcORJbtmzBmDFjsGDBAvzVX/0VWltbAQANDQ1obm4u+54vvvgiFi1ahJqamqKf9bJ161bMnj0bjY2NWLlyJQCgubkZ8+fPR0tLC+bNm4cbbrgBF110EdavX99Xn7a2tvB7AHj22WfR1NSEHTt2oKOjo+x6VZtisenl008/xX/+53+ioaEBgGMzmBSLx5e//GUAQJZl6Ozs7HuaUkk86uvrMXbsWADAiBEjMHz48OTP+qNiBFQvTqc6RsXi8vvf/x5f+cpXAAB/+Zd/ic2bNwNwXAaTYvEYPnx47knj1q1bMW/ePADAtddeiy1btvR9V26cSolHOXEaUjHKqsTrr7+ezZkzJ+vp6cl27tyZLVy4MPv1r3+dffe7382yLMtefPHF7NJLL82OHTuW/fKXv8zuvPPOLMuy7JVXXslWrlxZcK22trassbEx99/+/fsLfvfpp59mCxYsyLq7u7NLLrlkwM/688ADD2S/+MUvsizLsu7u7izLsmzhwoV9dbjmmmuyLMuyp59+OrvxxhuzLMuyxx9/PFu1alX4/VAlik0vzc3N2d/+7d/2HTs2g0MUj+7u7uyWW27Jzj333OzP/uzPsj/+8Y9ZllUWj17+8Ic/ZF/96lez48ePl/xZlukYZdmfRpyiuLz88st94+MHP/hBtnTp0izLHJfBInXeuv3227Pf/va3WZZl2Y9+9KNs48aNWZZl2a5du7Jbbrml73eVxqmUeJQSp6EUo6qaCi+++GIMGzYMs2bN6vu/stGjRwMAJkyYgK997WsYOXIkrrrqKjz88MO9C5Lcderq6tDS0hLeb82aNVi8eHHB/22qz/pz99134+GHH8a6deuwZMkSnHXWWTjvvPP6vr/wwgsBfPbop/+/d+/enfT9UKVYbHr553/+5wIty7EZPIrFo7m5GWPGjMHvfvc7/OY3v8H3v/99/NM//VNF8QCAQ4cO4dZbb8Xq1av7DEepn/XCMbr++uvx3//9338ycSoWl29961t4/fXX8fWvfx1//ud/jvr6egCVjRPAcSlGyrzVn4kTJ+LQoUMAgM7OTkyaNKnvu0riVEo8SonTUJvjqrog2L59O7Isw+7duzF16lTMmjULe/bsAQBcfvnlePzxxwF8pvOcffbZAIA9e/bg/PPPL7hOe3s7brrpptz1N2zYUBDgd999F9u2bcOaNWuwa9cuLF++HCNHjsx99thjj/WdU1tbiyeeeALHjx/HJZdcgjvuuAPXXXdd3/fDhg2T/+7tTNH3Q5VisQE+kwu2bt2KVatW9X3m2AweUTwmTpwI4LOFdK/uWUk8eg1xK1euxKxZs0r6rD8co+uvvx7/9m//9icTp2Jxqamp6euv//AP/4BrrrkGgOMymETjhLniiivw6KOP4rbbbkNzczPmzJnT9125cSolHqXGaajNcVVdENTW1mLBggX46KOP8Mwzz2DChAmoqanB0aNHMXnyZCxcuBBXXnklampq8OyzzwIA3njjDTz55JMF10ldtT3yyCN9/7700ksL/rgM9NmqVauwYcMGHD58GMuWLcOWLVtwzz33lFnjzw/FYjN69Gi8/vrrfbHpxbEZPIrFY968eXjhhRfQ2NiIY8eO4ac//SmAyuKxbt06bNmyBV1dXXjooYdw1113oaenJ+mzm2++ue86HCMAf1JxKhaXzs5OLFmyBMOHD8c111yDr33tawAcl8Ekmreuv/56bN++HTt37sSdd96JZcuWob6+Hg0NDZg5cybuu+++vmuVG6fUGN18883yt8XiNORiVBXhoQhNTU3Z6tWr5XcffPBBdu+99w52EYqydu3aU3r/U4ljM7QY6vEYiD/1ODkuQ4ti8RiIoRqnoRajYVk2BJ7TGWOMMeaU4tTFxhhjjPGCwBhjjDFeEBhjjDEGXhAYY4wxBl4QGGOMMQZeEBhjjDEGXhAYY4wxBiVkKuRdolQ+ev5szJgxBcdHjhwpOB41alS+QLQ7VE9PT9Fj4LO0u8Wuy+UAkNsv+vDhw7nfFCsXAJx22mkFx5y3musLFKabBIDGxsai902h/45eqlzAZ3tu94dj1X/Xw4Hgc3hHLd55LOW+3GZAvo04VUZ3d3fBsaovX4Pzn6vYcFmuuOKK3G9K4cUXXyw47s2x3p9x48YVvQbvp851B/J14/7/8ccf586J4qD2U+c25Xbna/buRd8fHos87tScwHXuzcRXCWvWrCk4VmXl+nA/5O8/+eST3DXOPPPMgmMeMynjju+rUsfwHMj96tixYwXHah7leHKbqLFa7dj0ZrAtdk+uC/d/7u9qjPFY5O2Cua8D+Xbn2PFYBfL9m/+ORHMdkI9Db1rzXtRcxufceeedud8o/ITAGGOMMV4QGGOMMcYLAmOMMcagBA8BayxKc2e9h3Uq1lyUbsd6COtaSkONtHt1Hz6HNTjWepQ+xJ/xNVN8B9WAtUulf/F9WSPmNpo6dWruGgcPHiw47r+Nq7oGUHzv8tRzODbcj5SGxvU9cOBAwbHylajrVAKPkbFjx+Z+w30oai/Vf1hT5WPVd6dNm1ZwzNq2ah+OA+uwrNWqscr9js/he6hzqgG3ifJE8fjl8c3tzOMByJe9HE9U5KEBgPHjxxcc85xQzrzDbaL6kSp/NVFjkscRl5O/V3Xn33Cbq3M4DjwfTpgwIXcOw+3FbarGHf/tZY+E+juT4k1R+AmBMcYYY7wgMMYYY4wXBMYYY4yBFwTGGGOMQQmmQjbHsNEByBvE2NjCRg5loGDYQKHuGxlblGGIDUFcNq6LMnuxEYXLwfcAdOKVSklJMsT1Of300wuO2YjICYWAfP3Y3KTMjNE5nGRD3ZvNPGz4KccMqJKVqMQylcCxVklHIpMRx031ZY43j1VldosMsWpMcX9mExVfU/X1KHmRSkQzGEZcHs+qv0eJh7i/q/7D9+E2U32X5zi+rupHXFZueza4qmvwnMBjVc0r1TYVcn9IuWdUBmXC5Lk7aj8g/7eI+2rK/MFtzIZAZaqN2kSZbpXRMAU/ITDGGGOMFwTGGGOM8YLAGGOMMSjBQ8Aau9LHWVNhjwBrjEq3Y22XtRCVEIk1lZSkQqzLcGIPvkY5CXSU9hltolQOKRprpBFzfJUexlp0e3t7wXFKMgwua2dnZ+43UUIr7jcqIQjrhnys+m+Kp6UUuH+oukZ9lb9XHgLWNlmnVHXlpCpcVqXlc/vwmOD+rzwlXD8uqxpn5SZZKQa3mUoaFZWD66v8AKxFc99VY/ejjz4qeg3VT3le4bbnsitPFLcJjxl1X7VhVyVwHJQWzv2Z2z1KXATk+zf/JiUZVpQwDcj3Z44/t6kaM1EiNnVf5ZtIwU8IjDHGGOMFgTHGGGO8IDDGGGMMSvAQsF5SW1ub+w3rktG7ner9UdZcWFNR78+yzsQai9JDo42WWINJed+WdSelXSmNqFJYl1W6JJeNNcR9+/YVHKv39KP6pbz7Ws475awZRhsApVxD6aFKa68Evl5K+3C7t7W1FRxPnjw5dw6PRR6r6t1m1u75WGmo0aY7KZvKsE7NY5XLAeT7ZjXg+Ks24nHF7czxVDk1WNvndlXvu3PZUvRgnld4bPKx8jtEuro6p9z33QeC66ruGXnG+BpqHua+yX05JYcMX1f5w/g3PHcpfxuzf//+gmOur5qrlRcnBT8hMMYYY4wXBMYYY4zxgsAYY4wxKMFDwLqN0sL5HUvWpFj7SMmnzZoj5wsA8tpNlA8eyOsu0Xu8SitjfYs1JqVLVvtddyCvS6l2jfYqYG1aXYPbMSWnOJct2ttAfcb9iNtV6XBcfo7NydjLgOuqvA9Rbv9p06YVHL///vu5a7Aeyscp7zaztq3agtuUxyLnWVB+gKj/q/uqMV8pHBulM3NZo/e///jHP+auwb4K7qvqvtzOfLx3797cOTxmorwiM2fOzF0j8p5U2y+giPKjAPlxdODAgaLX4LwOQH5M8JhRvrOOjo6CY55D1PzH8w7nbeBxp3wo9fX1BcccS3VfNX+n4CcExhhjjPGCwBhjjDFeEBhjjDEGXhAYY4wxBiWYCtlQw2Y/IN6o4Ywzzih6TUAb8SIicxvfF8ibCNkwoxKGMNFmLsrMdTKMOQquDxu+OKmUSirDZhY21SiDVLQBjjLvcFKQKJmRMq9FiaS4vkD5G4IMBBs3lTkuMtpGbQ7ky83XUGOK+z8b5FRcGDZzcfzVpivRdVWsy0lmFcHjQc1FbACM5hk1Z3A7c2xSNnNi45maQ/gzHkN8TZV4h9uZ+5rqR+XM18XgNlTtE226lDL/c7x5rk7Z/Izvo8zNKWOiPyohGN+Xr6nK6s2NjDHGGFM2XhAYY4wxxgsCY4wxxpTgIWDtS+l60QYbrNsq/Yk1JNZ/VOKSKBGL0nb4N1GSIQX/JtpASN2nGnD9lKbO7cYaU8rGU6z/8jU4MQ2Qjx/HSulffG/ue9yvVHwjfVfFV3kgKmHChAkFx6qukZafsulQ1KYpG6ikJISJ4hBp0KpsrBErvbfaCaOAfNmVJ4qTdUWbbLFnBMhr+ynJfqLxrOZevneUEEwliOL7crsr/bvavhueH1U5uc04dtw+al6Kkl2p+S/arE9tbhR5N7isPGcAcSI21UZqnkjBTwiMMcYY4wWBMcYYY7wgMMYYYwxK8BCkbJjCsPbBuo3S01m7ZF1G6bCs06RsMsQaEZct0jqBvMaYsnFRtTU3IF8/VdZIH4/eWwaAKVOmFByzTqXeMefNPDg2SkPlWLB2nbLxVPQuvmqjFN9IKaRosLyZCdelvb294Fjpx1FuBzVmuK78GzU2o02mUrR+HiMp752r8lcKzyuTJk3K/YbbINKQ1TzDMef4Kq2aNWL2A6g2qq2tLVoWbkMVKy5b9M48EOf7KBXuHyn5EjgO/D3PW0B+XHHfVnMBzzs8t6mxyeWPvBxqLuNzOC7Ku6D8ayn4CYExxhhjvCAwxhhjjBcExhhjjIEXBMYYY4xBCaZCNjsoU0pkGGMzjDKt8IYpfI4y7kXJPpSRL9rMgr9Xhjk21bEhUiWHSDFjlgongFFl5bbu6OgoOGZTlYpvW1tb0fuoup155pkFx5xEhM1/QJy8KSXJCv+G+4RKTFJtgxS3j0oQxGOG253LnZJwJKWuUfIbZebi+kRxUP0wMtWp2EdJZMohZTMzNmtxm3EyoHJMeCmbSDHKIMxljeZadd+oTyjz2mCjEvVwoiE28/GxmsvYzMv9XfV/HnvcL5XZk8cEtyH/7VJJlLjPcJtU06juJwTGGGOM8YLAGGOMMV4QGGOMMQYleAhYX0rZlIN1Oda+9u3bl7sGa12slyidLto0SellURIa1jJTtG4um9KUUnTGUmENqZzNYFgPVVoua3OsbSmNPGpH1Y+4n0QbgijvApcl2mgL0BvrVALfUyULYU8F67YpvhuuK7d5Shtzm6oxs3///oJj9q6kXIPvm+IRGowNwbis6h7RBkEpiazUvNEf1XejfpOymRcnKuL+nuJl4jlB1UVtglVNUrxqTMpGdVFcVJIh7ptcNtWm7Ang2PLcnTJWU3xEqvwp+AmBMcYYY7wgMMYYY4wXBMYYY4xBCR4CRumDrIew1sEajNKf+LqswSidl8+J3uUG8noPa0qs1ap3f7m+rNsorW8wPARR2YF8fZW/IboGx4+9GkqX499wm6j3fbms3PYcb3XfqKxKz6+2Vh31DwVvmMJjRHkf+N1m1o9VP1Sf9UeNb27TaLyrNuZ+xR4K5SEYjNwd3P/VPbj80Xv4KXkXuI+l9An2Wqk5kMcAx4L9Hgq+LvsSVBtVe7M2bi/V5tx3ub9HXicgP++wdydlnHFfVXNqSj6b/qi5LPIRqfqVO5f5CYExxhhjvCAwxhhjjBcExhhjjIEXBMYYY4xBCabCaCMfIG/CYcMJHyuTCpu/okQt6r6RYVD9ho0ZbHaaMmVKeI0Ug9RgELU7kC8r/4brr9qZDTFsXFEmMjbnTJ48ueBYJdko1XijTDVcfr6maiNOxFIp3B6qrpw0JTKHqTFTV1dXcMwJhJSpNjLApiTiYjMXj11lmONNtXgeUedU27gGaJMwEyW44bqkJCJjMxu3GZDvNxw/nmeAfH34unysTMNsxuZyqLJWe47jfqeSCvFnbDznZHbKIMjtxf1fGTc5djzvqDHD94lMpTyW1W+4HKqNItPwQPgJgTHGGGO8IDDGGGOMFwTGGGOMQQkeAtZHlAYXbaqSsukE6zJ8rDQr1r5SNoRh7ZJ1OT5H6WdRspso+U+1YJ1Z6UdcNtbU+RopWm60UQeQjxf3o/Hjx+fOYQ08SmiV4nfgvqb0/Gpv1MJaX+SNAPLlYj107969uXO4/6f0XY5vir+H/TzRZmaqjaPEXKqsrBFXg5TYRJtqpYz3aIyoPhf5RiZNmpQ7J9rQjGOn7httCKbmfJVYrBKiDeMUPEa4TKrcUVI1ldiH24fvq+bMSP8vZ2M6/puoxllKIiqFnxAYY4wxxgsCY4wxxnhBYIwxxhiU4CFgDUq9C6veQ+9PyoZB0bueKe8k8zukSpfkjSrY78D3Ue9yMymb/ajPKoXLqnR5/g3rnawpKl2Kz2HNWOUDYN2NNTLVj1j/4naN/B6qbKwJKu1TtVslcHso3Zb7YaTDK22f68L9X3lKuA157Co9lOsT5TLgOQOIN11J0WGrAce6nH7I5VL6L3sTUvJ9cH/mmKuxyXMpX4P7WUpf5PqrOTBlXiwFnqdScsgw3Mbq7xL7lLivqjkl8gOoWHJ9+D5cPxUX7nccf/V3NNqIayD8hMAYY4wxXhAYY4wxxgsCY4wxxqCCvQxUPgDWEPk3Ke/+Rpqb0qn5HNaYlJbJRDnklT7EWu2hQ4fCskb6VzmkvJfKsWF9kPVPlSOC9UKOr8oZzu9Ms7al2oj7CWuA6hyG311nzVjpbtV+353vwf0DyPddPidlDxGOJR8rPZH3lIg8BYpoPxB13+i9chVblVe+UjgWKfuqRPWtr68P7xu9Mw/kxx7fV+25wf0oOlbX4HmS66/mhGqPmRT/C/tQuE25rsqDwnsG8LhS9+X5j/uQ+psY5XJgj4G6L8c/JQ9Lub4bPyEwxhhjjBcExhhjjPGCwBhjjDHwgsAYY4wxKMFUyGYgZTpiMwebUvh7ZTpiUwkbJpT5J/qNug+XhQ0jXD9l5mJDDRuElMkkJbFSqXCbqfpGyX3YiKZMR9wHUpI3cWy43VMSj7B5i405qi+ySZKPlaG12glwuB4q+Q2b+6KNWVT/YTNYtCkPoM1L/VHmP45l1KYqLhxL7jOqrINhxOWyqf7OcCy47ypzF88BUXInIB+/lARBbD7lWPBcpdqU+wSXTZmGq52YSCWRi2BzX0qfijYVUufwvMr9Qc0p0SZ6XFbVntwm3M+UQbrcBHh+QmCMMcYYLwiMMcYY4wWBMcYYY1CCh4A1CZVQg7UM1s/27t1bcMxJawBg3759BcecdEIRJSZSumvkM+Dvo01ZFEoznjhxYsnXieBYKH0w2ryG66/0MNYYWQ9ViUu47VN8B9zW3I+4L6pNZbgskT4KpCWwKgW+p+rLXBduD9YPlS+F48/3VXoix5c9A0oPj8ZZyhxRW1tbtKzqvsq/Uyl8H5X8iPt3lHhNjZnIE5Lijzhw4EDBccpmXlFCK+Uh4evyGFFjtdqbtaVsGMXtHs1dqowchxQtPxqrKpYcB253NWcyXF/+u6J8F+WOGT8hMMYYY4wXBMYYY4zxgsAYY4wxKMFDwHqa0stY74neD1V6IeusrNuk3JfLqt6HZq2K9bNIYwLy7y7z+8Jqo5ZyvAgR3EZKY4y0fEZpUFx2vo/aVCnyDCjdjdst2uxIxYbjx2VXHoJy3oEuRpSXA8jXhTVFLrfqUxz/lP4f6ZKqf/B1uT7c/1V9+TO+ptJ7B2PMcL9L8cxwP+Pxr64RbRKVosFHvgsgfr+dx2aKRs5jSOn5vElQpXC/Uxo7zwfR5m5qPuT6q3f5mehvkcplMXXq1KL35b6txgzHm+ujvDrOQ2CMMcaYsvGCwBhjjDFeEBhjjDHGCwJjjDHGoAJToTJ78AYabDJis4cylLHZiw0VKtkPn5OySQ3/hs0dKWZGTmSTknQi2lSmHNicqYxMbPCMkvAoAxEbflLaiNuAk6xEpisgb/jhNlTJnri+kXlPXbdS+J7KdMlwm3KbK1Mht3uUuAXIG5FSzJ6RuYnHkOpjnHhs/PjxRcsB6DpXCtdPxT5KZpWyuRG3MxsEVTunxI9hM3aUaC1lrHLf4/kcAFpbW8OylUK0SROQrwvHJWUDKa4bt7kysipzYn9U321vby84jjYrU2buaFMlNT7K3UTPTwiMMcYY4wWBMcYYY7wgMMYYYwyAYVm1d6cwxhhjzOcOPyEwxhhjjBcExhhjjPGCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMfCCwBhjjDHwgsAYY4wx8ILAGGOMMQCGn+oCfBF47bXX8NZbb+G9997DqlWrMGLEiFNdJAPHZSjj2AxdHJuhSTXiUtYTgpaWFqxYsSL3eVNTE55//nl0dXVh9uzZGDduHHbs2AEA8jMA+PDDD7F8+fJyioG3334bDQ0NaGxsxOLFi3HixAn52UC/jejq6sLSpUvLKlt/rrrqKtx///2ora3F0aNHK77eQJQTl4HaxXGpLuXEZseOHZgzZw4aGxvxzW9+Ex9//DGAymLT09OD22+/HQ0NDbjyyiuxe/duAMCKFSvQ0NCApUuX4vjx432/H+jzgfiixKaXtWvXoq6uru/Ysake5cTl/fffR11dHebOnYu5c+eivb0dgOOSSlUlg6effhpLlizBmDFj0NTUhEWLFvV9pz4DgBkzZqC1tRWdnZ0l32/69Olobm7Gpk2bcM455+Cll16Snw3024hXX30VV199dcnlUjz11FOYP38+zjjjjKpcrxSKxWWgdnFcTg7FYnPuuefizTffxKZNm3D55Zdj48aNACqLzfbt23Hs2DFs3rwZDz74IH72s59h27ZtaG1txebNm3HBBRdg/fr1ADDg58X4osQG+OwPxfr16zFjxoy+zxybwSeKS2NjI1paWtDS0tK3WHNc0ih7QbBz504sXLgQl112Gd555x10dnbiyJEjGDlyJIYPH16wagYgP+uloaEBzc3NJZehvr4eY8eOBQCMGDECw4cPl58N9Nv+bN26FbNnz0ZjYyNWrlwJAGhubsb8+fPR0tKCefPm4YYbbsBFF12E9evX99W9ra0t/P7ZZ59FU1MTduzYgY6OjpLrWQqlxqVYuzgu1aXU2PR/5PfJJ5/gvPPO6zsuNzZf/vKXAQBZlqGzsxN1dXXYunUr5s2bBwC49tprsWXLFgAY8PNevsixAYAXX3wRixYtQk1N4TTq2FSPcuLy5ptvoqGhAT/84Q+RZVnf545LTNkeggMHDmDz5s3YtWsX7rvvPvz93/89Zs6cWda1zj77bGzdurXgs/b2dtx00025327YsAGTJk0q+OyDDz7Aq6++igceeKDoZ8U+/+Uvf4kHH3wQ3/rWt9DT0wMA2Lt3L6ZPn45du3YhyzK8/PLLeOaZZ/CP//iP+Nd//Vc88cQTeOmllzBr1qyi3//N3/wN7rjjjrLaplTKjYtqF8elupQTm//4j//A3/3d32HEiBH4wQ9+0Pd5ubGZMmUKampqcP755+PYsWN488038dxzz+FLX/oSAKC2trZvMuns7JSf9/JFjk13dzfWrVuHl156CY8++mjBd45N9Sg1LmeeeSb+53/+B2PHjsX3vvc9bNy4Ed/+9rcBOC4plL0guPjiizFs2DDMmjULbW1tAIDRo0eXda3+q7he6urq0NLSEp576NAh3HrrrVi9enXf/1Gpz4p9DgB33303Hn74Yaxbtw5LlizBWWedVfB/ZBdeeCGAzx5x9/93r54UfX+yKCcuA7WL41JdyonNN77xDWzbtg0//vGP8fOf/xz3338/gPJj09zcjDFjxuB3v/sdfvOb3+D73/8+GhsbcejQIQCfTWi9E+HEiRPl5718kWOzZs0aLF68OPd0AHBsqkmpcRk1ahRGjRoFAPjOd76DrVu39i0IHJeYshcE27dvR5Zl2L17N6ZOnYpZs2Zhz549ZV1rz549OP/88ws+S1m5dXd3Y+nSpVi5ciVmzZo14GfFPu+ltrYWTzzxBI4fP45LLrkEd9xxB6677rq+74cNGyb/3dvJou9PFqXGpVi7OC7VpdTYHDt2rG9yq62tLTAolRsb4LNJCwAmTJiAzs5OXHHFFXj00Udx2223obm5GXPmzAGAAT/v5Yscm3fffRfbtm3DmjVrsGvXLixfvhyPPfYYAMemmpQal66urj79/I033iiIg+MSU/aCoLa2FgsWLMBHH32EZ555BhMmTEBNTQ2OHj2K0aNH4/rrr8f27duxc+dO3HnnnVi2bJn8DPgscE8++WTB9VNWbuvWrcOWLVvQ1dWFhx56CHfddRd6enpyn918883ytzfffHPftVatWoUNGzbg8OHDWLZsGbZs2YJ77rmn3OY5ZZQal1GjRg3YLo5LdSk1NlOmTMFPfvIT1NTUoK6uDs8991zftcqNzbx58/DCCy+gsbERx44dw09/+lNcfPHFqK+vR0NDA2bOnIn77rsPAAb8vJcvcmweeeSRvnMvvfTSvsUA4NhUk1LjMm3aNDzwwAMYO3YszjrrLDz00EN913JcEsiqSFNTU7Z69eqSzvnggw+ye++9t5rFqApr16491UWoGo7L0MWxGbo4NkMTx2XwGJZlp+D5nDHGGGOGFE5dbIwxxhgvCIwxxhjjBYExxhhj4AWBMcYYY+AFgTHGGGPgBYExxhhj4AWBMcYYY1BCpsKnnnoq/E1tbW3B8ZEjR4p+37u/e3+6u7sLjk877bSC4xMnTuTO4XzinA9f0T/tI4DcvtUp9+X78DVGjhyZO4d381uyZElY1oi1a9cWLRcAHD58uOCY68Ptwcfquvwblded78O/Uft2cxvxfbldjx07lrtGhDqH79ubSbNcvve974W/4f7O25Z2dXUV/R7IjzNuc+7LwGe7J/Zn6tSpBccHDx7MnfPpp58WHHNOea7LhAkTctfYt29fwXFvauZeVB/isr7wwgu535TKXXfdFf6md/OZXrgfcrnUeOdrcB/jNgPy8UtJFcMxjsYmz1XqvtE4BPLlf/bZZ8OyFuPee+8tOFZ1578b3A+5n/L4UNdVsWO4fbjv8n2BfJvxvMP3VdcYP358wTH3KTWHclxWr16d+43CTwiMMcYY4wWBMcYYY7wgMMYYYwxK8BCMGzeu4FhpsPwZ61gHDhwo+j0AjB07tuCYtS6lMfJvuBxKQ2XthjUX1ocUfF/W3JROp8pfbbid1X1Zi2Z97PTTT89dg3Vlblelh7JWN2bMmIJj1sOAfLsp/0Z0jUgj5P6s7lspKV4W/s3+/fsLjrm9lF7I/Z37rtJQOb5cd74vkPczRHo4a+xAPL75moD2TVQKzz0pHiHuZ9EYUteI2lB9xn1XtWvkIWCdPcUjxP3qZGx7w/dMmVM4LqzDq3JzXXl+VJ4C7rscSzVnMtG8pPo6x4r9YOrvW7l/Z/yEwBhjjDFeEBhjjDHGCwJjjDHGwAsCY4wxxqAEUyGbFJSBghNGKMNMf9joApSXIIgNUGyiSkkYwrAxRSWMiJJZqPpF9y0HNpopY0pk1IvMT0C+7GyQUUY9NglxbJT5hfsN/4avoQxS3AbKjMcoc04lpNyT22zy5MkFx9zmyiDI5eZ+mWIG47GrTIUpRrz+KJNmZJhTZt5yEk+VSooJi9uxs7Oz4FiNd+4DfB9VNzbS8Vyr+mk0B3R0dOTOYXiOmzJlSsGx6kdq7FUCt4+6Pn926NChotdQceG68jhTfZfHamTMBWITNc+7nLhL3YfjoBL8pRgcFX5CYIwxxhgvCIwxxhjjBYExxhhjUIKHgJMhKP2ctR3WNljXUPoQ6z0pOh3DZVOaW7T5B5+jdBpOVMG/UZr6YCT3YD2MY6XKEmmbSlPl+HF9lR7KseCyqeQ93PZRcirVFzmJEmvTJyPJCvdVpf9zHLic3F5KY2edkhNTqfbhNk5pD64Px4GvofoDxzJKVKTKWg24HdVmTlHirUmTJhUcq4RBkWas6saf8TkpmjH3G95oSpWVrxElmgN04qBK4DGifGhcDt40j79X/i/+TUqSNW4zvm5KP+Xf8FhV9eU4sHeHNyYD0jwjCj8hMMYYY4wXBMYYY4zxgsAYY4wxKMFDEG0GAeT1D9YcI30UyOtYEydOLDhWWh+XLcVnwHoQ64EpfgCGtVzVRtXW3NR91D0izZh1KaU7c5ux3qt0aP4Na3fqHP6M65Oiw7KuyDqc0hUHe+Mppf9z+/Dx+PHjC45VfogoD4Fqnyguish3w6icG+wZKCdHQjXguUflVOB25PrzXKT8MNHYVOfw/MU+k3Le/Y/eu1f3Yd+B8oREuWZKheOg4sLl4jhEngIg3phL/Q3hvsrjSnmEor+b/PdOxYU/47KruHAOiVT8hMAYY4wxXhAYY4wxxgsCY4wxxqAEDwGjdJkoHzprVOpdWH5/MtLxgDgPv7oP5/p+//33C45Zp0rxTHB9lfapdORKid7bB/L6FpeDdSilVXPbRzockM8jEelwQL4dWf9MyQcf9T2lEUZ7U5QK31PprVwX1gu5/fgYAPbv319wzHVTuiTHKkXL57LyGOI+w74cBbdRig5fDTgWKft/RB6pPXv25K6h9oToj9LI2WeRsp8Lz3E89/A1Zs6cmbsG+ypYm1d9pNrzGbexuj73s2isq/0BuC7c7/bu3Zs7h/08HAflS+J9Frh+HKcZM2bkrhHlVVC+Ou9lYIwxxpiy8YLAGGOMMV4QGGOMMcYLAmOMMcagBFMhGyiUaSFK3sNmEGV2KmcTmsh0poxrbDxh2tvbC46nT5+e+w2XjTeq4KRKgDbAVUpkdgJisxq3kWozNh2ltHM5G9Ow4Y3rx31RGdHYZMV9Uxmzqr2JDpvS1D257NFmP8pkxX2ZjUxqXLIRi9tLjU2GzYxcdnUN7jMpSaYGw1SYYjJV8Sp2jhrvkdFaGXGjTYVSNhlSY6I/KokOf8b1UbFR16kmKWOG24vHiGpjTl7Ef0PU3zfuz5EhGIgN3ymG0Wh+Z7MjoI30KfgJgTHGGGO8IDDGGGOMFwTGGGOMQQkeAtY6lEYRbXSRsikD6zKcLEYl3WH9J2WzIy4/JyJiD4HSofg+rDsprS/FE1EqrEsq/ZDbjduZ21AlVeJz+FglM2JdmeOp2pXLwroa9yPV77idWd9TurrqW5XA/S5lkxXWtqMEUuoa3KYqYQpfh+OkYsn9O0VDjcrKZVNtpBKvVAqXQ81n3Fd5XKUkd+L7cF3UJjRclpRkPezv4b7HYyjF75Pis1B9qxK4H6r5IUpuxu2l6sr9P9pQDcj3TZ4vVP/n9uE2jPwPqmwcF+WHKzcufkJgjDHGGC8IjDHGGOMFgTHGGGNQgocg5Z1z1lBYH2H9hN/bB2LNRb2THGl7SndlXYbfp015L591OdYclS6ZorOWCpdDaZn8GevBrLupzT04FnxN1v+A+F38lBwRvPEM51BQG9Ow3hfljFDnVAr3S/WOMff36L1k9X409++Uza74Otw+qk35upEOqzTn1tbWgmPW0NX4UNepFPbIsIcIyLdb5EtRWi7/huuizom0fPXuf+TP4rlJzaORv4fzWwBx7plS4bqq+YHHaZTbQsF1S9mIjD/j9lLzLseB7xtt5AbEuUpSchek4icExhhjjPGCwBhjjDFeEBhjjDEGXhAYY4wxBiWYCtmEk2I6igw06hpshuCEGykbObBhSpmqIkMUmz2UMZHNHnwfZW4ZDFMhG1FUWTk2nACjo6Oj4FgZ7KJkF8pUk2KsY7hf8DlsZFL35cQ63CeUGYr7WqVweykDGW8gw7/h/s6GSiA/zlI2uuE242M1ZrgPRWYudY0oYU61k0MNREqCLK4Pj10eI2ps83W5b6ck5kpJtBMl1uGyqT7B53DfU8nKqj2fcZ9S/YHvySY7bh+OtbpGSrIznkN4rKpYRmZO7v9qjuV+Fm32BpQfFz8hMMYYY4wXBMYYY4zxgsAYY4wxKMFDwPpJyuZG0eYg+/bty10j0imVzst6OOswSpdhvYfvy5qc0pQ4QQi3UcrmTdWA21XVt9SEN6rsEcp3EGl1SuuKksbs37+/4PhLX/pSeN8oqQygtdlK4FizL0V9xnXjOCgfAvdljqXSNiOdVSW/4XP4uuzbqKury12D2501da4/MDiJiTg2KfpvtEGYguPLY0RtZsO/4XKoOTDyO6Qkq+IxwnO88h2oflIJXPdy/BJ8jhrrXH/uY+rvW5TMTiVmi9qUx/fkyZNz14iSDKm4lDuX+QmBMcYYY7wgMMYYY4wXBMYYY4xBCR4C1jqUXsKw3snvfqv3QyP9P2VzF9aU1DvnrP/xMWtZKRv3cFmVZqw2CKkU9gMozYm1LNYuuQ3VNaLcDSpHBLcJ63mqjfg3rE2zp0C9Hx1tbpTyvm+lpGjM7e3tBcfcv1mnVG3MGiL3MaXBR+9QKw2V76P8DP1R9ee+ytdI6Q/VgNtE9Qdue64Pl13VN9qsTXkmOMYpcwa3a1RWzn+hzuGcF6ofpWwkVArcD1WOEe4jrJcfPHgwvE+0QZKC50Q+R/lBIp8Bx0F5O6LN+1TZlScmBT8hMMYYY4wXBMYYY4zxgsAYY4wx8ILAGGOMMSjBVMiGiZRNOdjswEYlZRhhEwYbD1M2x2GTSUqCIDYupSQ3Skk8wyiDY6VwuyoTFt+XzUwphjqORUo7R5uqqPtGBik23iijHRuL2BCk4jDYG+ukJFmJ2kcZpthUyXFSpsvIdDRp0qTcZxwXnhM41qofsqkqJdlPtKlWOZSzQQz37yjpEJAvOxtkVZ/gNok2twHy7cZlZSNuyoZgnDRKlbVc89pARBt1Afmyc//muKQkpksx77K5s7a2NrwPj0UuO8dJzUtc1sjMrs5JxU8IjDHGGOMFgTHGGGO8IDDGGGMMSvAQRAk2gLwGVc7GHqzdsD6UklCGf6OSZ7DGwuewTpeyYRAnmVH1U59VSsrmL6w7sVbFGqrS07mduQ1VO0eJiTo6OnLnsHbJseCyKe8CX4M1QqUZp3hASoH7rkpuxdpvlGRKaYMcf05Ck+IP4TipDXSiZF3sb1Bl5ThwP1TadopvqFQ4Nim6PI/dyNsCxL4b5Zng+LH+zwmD1HU4NhxPNVZ5Y522traCY7VZVbXns8hzAeS1e45d5EMD8u3F/VD5bvg+HEs1l6lEesVImXe5jVT9yvWq+QmBMcYYY7wgMMYYY4wXBMYYY4xBCR4C1pxYpwHyehJrX6xjqQ1UWB+LNuUB8loOay5Kh+L68DVY/0x5r5O13GpuOlEKSg+NNvPg9lB6OseP20zdl8/hcij9K2oj/l55KLgfsaam3qmONuspFfYxKC2cxwi/6xxt3KWuy7FVcWGNOWVDJL5OlP9CtTHHm9tceQhUH6k2aqxyX+U2iTxTQF4T5jZRfY7bmbVp1tCBfD9hPwf7A5TvhtuZ76Pm62rHhuuu6spx4TblOKXkF0nJQ3DgwIGCY25Dlbsj+rvCXgU1ZqKcOWqslrshmJ8QGGOMMcYLAmOMMcZ4QWCMMcYYlOAhYE1CaZmsf/D7kqyxqXclWXNmrZN1HFUW1lzU+7OsfbE3gfVB9S43lz/lnJORlz3FQ8DHrKEpDYq1Kr6P0mEnTpwoSvz/SdFduazcR5RGyJ4I1mpVHFjzrxTW9pVezPXncnMcVHtyXKIxBOTHDGumShtmvZPHFfd/1R/4s5RcFoOB0mqZKP8B14V1eiAfPx4zU6ZMyZ0TaeQpeRlU2/dHvWfP1+V4DsY+LAyPZTUPcTn5XX8e2+zTAeK9DNSYYXg8q7mM/65w2XnOSfn7wPdR/UH9DUjBTwiMMcYY4wWBMcYYY7wgMMYYYwy8IDDGGGMMSjAVsrFFJYyINmZRCYKia6hEJVHZ2CDFSVeAvHmDf8NmH2X2UNftjzKmqIROlcLGG2W8ZCNKZNxTRkyuL9dFGbW4D3A5VHz5OlxW7nvKvMnXZfOO2oik2kmjojIA+eQmkXlXGYjYAJhiZoyMauocNm9Fm8goYxNfg+OgEo+pjZYqhWOjxjeXlduEY6X6Ml+X21mZN6OkMuo+0QY4fKwS77DRkMdZSlKsSuExwhvGAXFc2FCZkiCOr6k2kKqvry845vZSCYK43aPEbCr20QZnaqyWa871EwJjjDHGeEFgjDHGGC8IjDHGGIMSPASswyj9nH/D+hJr6kpPZ+2DNRa1KQfr3f/7v/9bcKw2neB7s/7D3yuNMUq8pLS+amtu6ppKy+KysIbM2u2+ffty1+D6pWx+wppYyiY6XFbW91hTVjos626sRSp/R4pfpRRSkoOwrybq/wruh+x/UX2O25DLqtqH48ux4/sobwd7BrjsSqeOkuyUA18zJVEPJ7iJYgfk68PjQc2BPF9xu/N4APLzIvtuOKmQ0sjZVxDFGyg/Ac5AcD2U74bbjGMZJVAD8rHluCiPHLcZjxE1f/Bnkf6v/s7wffg3yruVsqGTwk8IjDHGGOMFgTHGGGO8IDDGGGMMKtjcSOkW0YYykacAyOsjrIUoTZV1J9bgUja/YZ2Jy670IT6Hf5NyzmCQstkF627chkrL4nNYY1R14zZgbU7dh3Vm1gQj7RrI+x1Y/63mhiADwe2l3qfn/sx9l3VLpePyfbjvqnHG/T9F/4+0/Kh/KHheUfdQZamUaNMhIO+r4VhE+TKAOD+Aui+PRY6f6kece4R9VXzNlM2ruGxqzFQ7dwf3S5UvhNuDY8nto/5WqXwX/VFeNb4O3zcldwfD56hypfwtYsrdiMpPCIwxxhjjBYExxhhjvCAwxhhjDLwgMMYYYwxKMBWyUUOZHyKzF5tulCEl2jBEmY7Y7LF///6CY07UAuTro0wk/Zk8eXLuM07ew/VX1zwZpkJlGOJ2ZMMXm1CUeY1jU87GHHxfZTxj0xD3G/5e9aPI0KqMRtGmMqUStRcQb17EJivVXpGZV5nyuGwppsLIiMh1UXNElBBMxUD1xUrhuqg+xGMiMkCq+ra1tRW9r4LnDY65MvdxWaPNq1Q7c334mqrvVTtpVMrfGSYyHSsTJps7+ZyUBFmRyRCIjZpsdj548GDuGpF5VfUplWgrBT8hMMYYY4wXBMYYY4zxgsAYY4wxAIZl1d7RxRhjjDGfO/yEwBhjjDFeEBhjjDHGCwJjjDHGwAsCY4wxxsALAmOMMcbACwJjjDHGwAsCY4wxxsALAmOMMcbACwJjjDHGAPg/o3eKtIR5Qp8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 15 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "im = ivim.io.base.read_im(im_file)\n",
    "\n",
    "fig,axes = plt.subplots(n//5,5)\n",
    "for i,ax in enumerate(axes.flatten()):\n",
    "    ax.imshow(im[...,0,i], cmap='gray', vmin=0, vmax=1.2)\n",
    "    ax.axis('off')\n",
    "    ax.set_title(f'b({i+1}) = {int(b[i]):d} s/mm$^2$', fontsize = 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing\n",
    "Some degree of preprocessing is usually required before model fitting. One example is directional averaging, meaning that data acquired with the same b-value but different diffusion encoding directions are averaged. While the simulation has not explicitly used different diffusion encoding directions, there are repetitions of b-values and the data can therefore be used to show this preprocessing step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "outbase_avg = os.path.join(folder,'comb')\n",
    "ivim.preproc.base.average(im_file, bval_file, outbase_avg)\n",
    "im_avg_file = outbase_avg + '.nii.gz'\n",
    "bval_avg_file = outbase_avg + '.bval'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAACSCAYAAAA6otLhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAFmlJREFUeJzt3WuIFuX/x/GPZuaaZgVG0YGCsMODSsqMZNvoYGVkPeiIdCCC6ABhGPigfBIRFZ0gCDtTkhCSQUK/jSgrMgpKKaksQioCM4rFdVfXXOf/oP9u7t6f7+713bl31+z9gj/8/rMz98xc1zVzX9595jsTqqqqBAAA/tMmjvcBAACA8ceEAAAAMCEAAABMCAAAgJgQAAAAMSEAAABiQgAAAMSEAAAAiAkBAAAQEwIAACBp0ngfwH/N+++/r88//1zff/+9li9frgMPPHC8DwnjjDGBPowFSOM3Dsb9F4K1a9dqyZIlDcvXrFmjV199VZ2dnZo7d66mTZumjRs3SpJ++eUXLV68uNZ+lyxZotbWVi1atEi7du0adv3Ozk4tWrSo1j4l6YILLtDSpUs1Y8YM7dy5s/bn/ZsN1/dffPGFWltb1dbWpmuvvVZ//fVXU/peklauXKmZM2cOuyzSrPEg/ffGxHD9vnHjRs2bN09tbW26/PLLtX379tr97sZS6bIS3B/yhhsHffa+LuuOgz179ujmm29Wa2urzjvvPP3444+S/PdB9jtC+vePg3GfEEReeOEFXX/99WppadGaNWt09dVX9//t2GOP1ZYtW9TR0TGiz16/fr22bNmijz/+WKeeeqpWrVo17DbvvfeeLrzwwhHtb7Dnn39el1xyiaZPn96Uz9vf9PX90Ucfrfb2dn344Yc68cQT9dZbb9Xue+nvm8KqVat07LHHDrlsKM0cDxJjQvqn30866SR98skn+vDDD3X22Wdr9erVtfvdjaXSZSW4PzRP3ziQGq/LuuNgw4YN6unp0ccff6wHHnhAzzzzjP0+GMl3hPTvHwf7xIRg06ZNWrhwoebMmaOvvvpKHR0d2rFjhyZPnqxJkybZf7W1traqvb19RPv79NNPNX/+fEnSpZdeqnXr1jX8fe7cuWpra9OyZcskSe3t7brkkku0du1azZ8/X1deeaVOP/10rVq1qv/Yt27dOuzfX3rpJa1Zs0YbN27Un3/+OaLj358M1fdHHnmkpk6dKkk68MADNWnS3/+Fq07fS9Lrr7+uq6++WhMnThxyWZ/RHA+S/pNjYqh+3/vn0e7ubp188smS6vW7G0uly/bmxoLE/WGkhhoHkr8u64yDY445RpJUVZU6Ojo0c+ZM+30wku8IaT8YB9U4++CDD6p58+ZVe/bsqTZt2lQtXLiw+uyzz6rbbrttwHo333xz9fXXX/f//++88061bNmyAets3bq1amtra/i/P/74Y8B6Dz30ULV69eqqqqrqhx9+qG644YYBf7///vurt99+u6qqqurt7a2qqqoWLlzYf7wXXXRRVVVV9cILL1RXXXVVVVVV9dRTT1XLly8f9u/4R2nf//TTT9W5555b7dq1q6qqen2/e/fu6oorrqh6e3urM888M1y2N8ZDc5X0+7vvvludccYZ1Zw5c/r7sE6/9xk8ljLLqsqPhapiPIzEcOMgui7rjIPe3t7qhhtuqE466aTq+OOPr3799Vf7fTCS74iq+vePg30iVDh79mxNmDBBs2bN6v9X05QpU4bcpqqqhmUzZ87U2rVrh93fYYcdpm3btkmSOjo6dPjhhw/4+1133aWHH35Yb7zxhq6//nqdcMIJ/f9KkaTTTjtN0t8/Q+79v/v+e9Rwf8c/huv7bdu26cYbb9TLL7/c/y/HOn2/YsUKXXvttQP+xeGW7Y3x0HzD9fvFF1+s9evX69FHH9Vzzz2npUuX1up3yY+l0mV9Bo+FBQsW6Ntvv2U8jNBQ4yC6LuuMg/b2drW0tOi7777Tl19+qXvvvVdtbW0N3wfZ74j9ZRzsExOCDRs2qKoq/fjjjzriiCM0a9Ysbd68echtNm/erFNOOWXAst9//13XXHNNw7pvvvnmgA4955xz9Pjjj+umm25Se3u75s2bN2D9GTNm6Omnn9auXbt05pln6tZbb9Vll13W//cJEybY/903UIf7O/4xVN/39vZq0aJFWrZsmWbNmtW/TZ2+/+abb7R+/XqtWLFCP/zwgxYvXqzJkyc3LHvyySf7t2E8NN9Q/d7T06ODDjpI0t9t3xfoqtPvbiyVLtvb4LGwYMEC/e9//2M8jNBQ48Bdq08++WStcSD9/Q9CSTr00EPV0dFhvw9OPvnk1HfE/jIO9okJwYwZM3TFFVfot99+04svvqhDDz1UEydO1M6dOzVlyhQtWLBAGzZs0KZNm3T77bfrlltu0UcffaRnn312wOeUzhJnz56tI488Uq2trTruuON03333Dfj78uXL9eabb6qrq0u33HKL1q1bp7vvvruZp4z/N1Tfr169WuvWrVNnZ6cefPBB3XHHHbruuutq9f0jjzzS/7/POuusAV/80TLGQ/MN1e/vvfeeHnvsMU2cOFEzZ87UK6+8Ikm1+v2NN95oGEt79uwpWnbdddf1f87gsSCJ8VDDUOMgulbrjIP58+frtddeU1tbm3p6evTEE0/Y74O+DFPpd4S0n4yDsf+vFGXWrFlTvfzyy/ZvP//8c3XPPfeM2bGsXLlyzPaFfavvHcbD6NjX+z3CeGguxsH4mVBV/8HfqQAAwAD7xGOHAABgfDEhAAAATAgAAAATAgAAICYEAABATAgAAICYEAAAACUqFX777beNG0/ym+9dhrHP9u3b7bp9JUr35t49HdWZ3717d8Oy6HWR7r3Sbv9RaYaurq6GZX1vRhvMnYP73Ghf7nyjdV0/nHrqqXbdOlwlsKhf9uzZU7yuO68DDjigaJkk7dixo2FZNDYdd1yZ8+rt7S3el7s2on255X/99Zdd17Xh+eefX3xcGXu/q76PaxcpN44d17aD3y/Qx7VttC+33I2ZqL3d+Wbuh+6+lbk+Im7dvkp6zXTnnXcWr+vuhS0tLXZddy276z7qF7du1H6lx+XWk3y/RtxxlY6LaPuIa5u+ap/D4RcCAADAhAAAADAhAAAAYkIAAACUCBW6oEIUgJgyZUrDssmTJ9t1S8MSBx98sF3e3d3dsCwKOLnlLnASbe8CiFG4xYUNXTglCiK5tu3p6bHruvYeDe6Yov7LBLRKQ4VReM+FzKJ9ueNy+4+CRJmwYmkYKgo9ZYJnmdBRXZlrxh1vdC9wbeNCXtF9J1peelzu+ora1S2Pxow736gfnUwwM+qHsRBdn66tXHhQ8tfXaAUw3T3abe/u+9FxZcaLE91fXCA+CtdG11cJfiEAAABMCAAAABMCAAAgJgQAAEBMCAAAgBJPGbjkYpTqzaSQSxPmUWrcJS2j43JpUZcsjlKlTpQedWWOXRu69GgkSo+6Jy1Gg2vr6PwzpYtLk/vR9m5sRAlcJ1OGNnNcpSWRo+1dajsql5pJ2NfljiFqL9e2madNonUd146Z4yp9CimSSaO7z43GbObJgczTC3W4sZkprZ4pJ5x54qj0KYWhlg/WjPLsbhyXnquUK7FfZwzwCwEAAGBCAAAAmBAAAAAxIQAAAEqECl0AIwp2ZAJOLmzhQhFR2MIdQ6aMqwuhRCVIXYgjE/Jy60ZBIrevTJBmNLhjjdoq815yF85x55oJXUVBLHdcbrxF27vzzfRhhjvWTFnu0eL2FV0Hrm0zAWH3uZlrLrpmonE7WBTkzYQlXRu4c8iELaPjis632dzYjALSmbCnW54J17p+yZQTztx3Xcn4aFxlvtMcNzYyJcBL8QsBAABgQgAAAJgQAAAAMSEAAABiQgAAAJR4ysCJSiS6BGiUgi5NX2bSwlHq26UvM0lP97lRkrylpaVhmUtBR09kZPY1VmVr3X6itLBbt26/RGPIjbdMOWEn2pc738wTME6UmHap7ShZnBnHdbk2jMaga6/M0yJ1n6DJPAWTeTLG9W3dctfNeFKk7pMtpVy7ZsrT79ixw67rkvuZEsFO5imgumX3M+XZ3TUTXceZ74NM2fbB+IUAAAAwIQAAAEwIAACAmBAAAADVDBVGIZrMu75doCoTjHHhvah0Y2kZ1EwgLROOyZQudutGwa0oaNZsmSBPpoSqawMX5InOP1O+ujQ41tPTY7d3Yz4KBZaWXI1CS6XtMtQxjAbXXtG9wI2ZaGy4vvnjjz8alkXn6pZnAp9u+yiU6M43c31Mnz69YZl7573k72fjXcbccfc8ybdhdJzuXH/77beGZVG/ZkrRu/bOlC52xxCNTXcMRx11VPG+XMA4WpdQIQAAqIUJAQAAYEIAAACYEAAAACVChS5sEYVgXCAr865wJwqhdHd3F20vlYfvMhX1MqFAdw7RO8Td50ZtFX1Gs2UqhmUqDZYGeaJ9ue2jgJMbmy7wE41td1yZqoZOtH3m3e6ZEOdYcucQna8Ljbp+iAKfThQoKw1hZsKaUR+45ZkgbGkYWhq7SoXuXtiMSouuDdwYyFwzmaBl5vp2fRCNAXfvisam49ogGi+Zzx2MXwgAAAATAgAAwIQAAACICQEAABATAgAAoJqli7u6uuxyl5CPkp6ladso6ZkpF+o+1yWWM0nd6CkDl/x3606bNs1u70oij3fp4kziOlPe1vVtJinr2jUqKe3GizuuTLI46pfS7TOleKN165QrzXLXUSZdnfnczs7OhmV1U+OR6MkSx/V5dK6uz9xTMJmnFFzJdmnsnjjKHH/d+3nm6QX39EPUr64PM9ecEz3d5I4rM47d9tH3VJ0xwC8EAACACQEAAGBCAAAAxIQAAAAoESrMvGvchTUyYYvMei5cE4W8XIDQhbGiQJtbNxP0c20YBTOnTp1atP1Qx9BsLsQSHVPmXeGOCwJFIRzXX5nQ0miFJd3nuv1nSjJHQaJmlI0tlSlh7dogE4B0nxuN90zbOu7+0IzApzve0kCblAvNjlUJ60wZdzdeMtecE10H7nOj+7m7R2fGthOV1Xb3CPe5USDQ3Q+j88qEYxuOacRbAgCA/QYTAgAAwIQAAAAwIQAAAGJCAAAAlHjKoLTEcLQ8Sk+6dV2Cdfr06XZ7l7aNUqHuHDLJeZfqjPZVmo7OlJ+MShRnSi03W7Rvl6qN0uEunezSulGyujTNL41O4jnzlIdrl2Yk2TNJ6Lpcn2euuUwa360bXQduX9ETIG5fbhxmyvFGT0K54y29F0XbZ57qGCtRwj5TWr103eg8S9P8Qy0fLBoDdcscZ55oyIy3Ok+d8QsBAABgQgAAAJgQAAAAMSEAAABKhApdqKEZZVTdZ7jt3XvRJR9CyYRzuru7G5ZFQaQoROm4ksSl5YwlHy6JSlVm3qtdhzvWKHDj+iDzvvdMiDWzbun71qMyrG68ZPqwbhgtU7Z3tLhzyIyDqG9Kx3F0f8ncdzLhTscFXKOSsa7PMqHbTEhsrEKFrg+jMejONQoIu8/NBIFdifxM+5XeHyR/XtE92oXE3b4y9/LoOqozBviFAAAAMCEAAABMCAAAgJgQAAAAJUKFTibIlHk/vQtWRMEr97mZ8J1b1tLSYrcvfa+5VP7O90wVragNo7ZpNhcQiyqLZd73XhqCicKimZCba1fXflFbu2ON1i0NkUZBIne+mffAjxbXtpnKkJnqf5kga2mIM/qMzPh2MuPAydxPo/OqU6UuI3N/cvfjqE1c+K70XhodQ/R94I7B3QsylWujNnDnkKkw6/YVjYE6lWv5hQAAADAhAAAATAgAAICYEAAAADEhAAAASjxlMHXq1IZl0fuvXQI1Kj3sEv1u+yhp6pZH65a+fzoqq+lSqVFatrQcbtSGLi2cSWePBrf/qK0ziWmXis2k5t32UdratVXp++qj44rGgEssu+so2pcbQ1Ffj9WTJlL95H/mPe6Zd8ZHx1BHdKyZfbmEuGvDaMy68ZUpEzwaMveczFM87n7i2jpzfUdt5e69rq0zTwFF48KdlxsD0Xhz55ApBV+KXwgAAAATAgAAwIQAAACICQEAAFAiVOgCGFF4wQUFo3LALmzhQhxR+C7zrmx3vG77adOm2e23b9/esMy9fzv6XBcYyYSDxjtUmAmIuT7IhHucKNzj2iUK+rl13bLMu8brBuqi8eq2j8qSjlWYTMqFS12gLnOsrm3qliOWyu8bmeszU57dHWvULm7MZMbcaHAlhiOuv6J+KS27G13fmbCm25c71q6uLru966/ouNw9zh1X1H9u3ei7J1PquWH/I94SAADsN5gQAAAAJgQAAIAJAQAAEBMCAACgxFMGrpxvlPZ1ZYqjJLVLIbsEbZS+dEnRbdu22XVLyyS7Y5L8+WZKxrpkbuYphUy6el8UpbBLn5KI+sX1YbQvN15cWjhKJmcS8u4Y3BiM+i9Tvvqggw4qPq666j79ELWhO19338iUTI/uG6VPGWSeosmclxvzUXLffW7mvEaDS9NHY8Cda3ScpU/8ZEr8Zp4Y6u7ubljmyo1HontUaSn26Fhde43GU2f8QgAAAJgQAAAAJgQAAEBMCAAAgBKhwqgko5MJNZQGCKPwnQswHnLIIXZdF0RxZR4z72uPAielIZCorUrf1R197mioWyI4E+5x7ZcZg3VLuEYhPRecypTKdttnAoyZ0qajJVO+OVMO2Jk+fXrxvtw9IgoglgbVomPNhArd57pl0ZgrLbk+lkpLgEfLo7Zy7XrYYYcVrSfl2rU0tBuVAi4tfRx9ruvDTBi6NKiYwS8EAACACQEAAGBCAAAAxIQAAAAoESp0wYwoyOTWjSo4ubCEC1ZkqgdmghlOphqdq+AolQc+mhEOyhxvHZlQnzuvTCUzJ1OtsrSvI1G/uIBRFFxz67rPzVS4y4SxRovrhyiI6443ai/XNpnKc9E9wnFtm6m+55ZH51UaDs1UnhvvUKE7puj4M9dMqcy9ILpHu6qhpZ8p5drAHW/mvu2u7+i4MtfBYPxCAAAAmBAAAAAmBAAAQEwIAACAmBAAAAAlnjJwSc2opKN7p3S0rktPuvRmV1eX3d490ZApjZpJw2fet16aCo1SqS4FHSWLo7LOzVaaApd8f0fruoR6psRvZl/uHNz+M8n/KOFfOraiMeCOIVPmeLTULUsdnYNL+ZeWf5Z8G9R9+iFTwjpz38mU8y0t7y6N7dMmg0XH5O7dUbuWfh8044kMNw7dOUTH6sZL5umm0lL60edG69Yp284vBAAAgAkBAABgQgAAAMSEAAAAKBEqdEGFqERiJjBTun0U7HDlJzMBRHcOUVjDhUCiMJVbvnPnzqJjksrDb2MpUyI4Cv04peVpo7BM5n3xrg1dkCkKZ2XawB1vphyxO6/Me+BHS6a9nEy41H1upr0y952DDz64YVmm5Hq0r9JwaTNK5Eb3rmbLfB+4sGg0BkrL7kb3XRf0y4Tv3OdGZZbdeI3u56Vh5Oj8M4FuQoUAAKAWJgQAAIAJAQAAYEIAAADEhAAAACjxlIFLakapWJf837Ztm133kEMOKVo3Sop2dnY2LIvSvqUpf1emWfIlgqPjcvsqTZ1H606dOtWuG33GWMiU7R2tFLY7hiiFXfqkSaZEcOaJBJdizpTizbT3aKlbujg639IkdaZvMml8lybPPD0RHZf7DHePjO477ljH+2mTzBMdmeN391jXrlFbufESXRulT3JFbZq5lkufFImONfP0SKZ88mD8QgAAAJgQAAAAJgQAAEBMCAAAgKQJVZ0EAgAA2C/wCwEAAGBCAAAAmBAAAAAxIQAAAGJCAAAAxIQAAACICQEAABATAgAAICYEAABA0v8BxJUUzdPyxpkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "im_avg = ivim.io.base.read_im(im_avg_file)\n",
    "b_avg = ivim.io.base.read_bval(bval_avg_file)\n",
    "\n",
    "fig,axes = plt.subplots(1,4)\n",
    "for i,ax in enumerate(axes.flatten()):\n",
    "    ax.imshow(im_avg[...,0,i], cmap='gray', vmin=0, vmax=1.2)\n",
    "    ax.axis('off')\n",
    "    ax.set_title(f'b({i+1}) = {int(b_avg[i]):d} s/mm$^2$', fontsize = 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model fitting\n",
    "The final step is to generate IVIM parameter maps. We here use a segmented algorithm with a b-value threshold of 200 s/mm<sup>2</sup>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "outbase_fit = os.path.join(folder,'fit')\n",
    "bthr = 200\n",
    "ivim.fit.seg(im_avg_file, bval_avg_file, regime, bthr, outbase=outbase_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the resulting parameter maps for a single slice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAGDCAYAAACYxjKAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAMDBJREFUeJzt3Xl0ldW5x/HfYQxDBCIBiYQAgZRBFBRxAAcUahWF9qp0Wa1gSq91vNCqdQJFFKeKQ51FBKp4xQJVqYJA44DihAMiMmgYDTWKYAIEgeTcP1zJJc27nzfnsJPjCd/PWq4F7+O798457znn4eR99hOJRqNRAQAAeFIv0QsAAAB1C8kFAADwiuQCAAB4RXIBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVycVPxNSpUxWJRCr+S0lJ0SGHHKKBAwfq9ttvV2FhYaKXCFTLc889p549e6pJkyaKRCL6+OOPE70kQFL877NTp07VzTffXLuLTXIkFz8xTz31lJYsWaIFCxbooYceUu/evXXnnXeqe/fuWrhwYaKXB5i++eYb/fa3v1V2drbmzZunJUuWKCcnJ9HLAiqpzvvs22+/raefflr/2SHj5Zdf1osvvpiIZSeVCL1FfhqmTp2qiy66SO+//7769u1bKbZhwwYNGDBA27Zt05o1a9S2bdsErRKwvfXWWxowYICee+45DR8+PNHLASqJ5X12z549uv322/XJJ5+oV69e2rZtm7Zt26bWrVtr3Lhx6tq1a4J+iuTANxdJoEOHDrrnnntUXFysxx57LNHLAQKNHDlSAwYMkCT9+te/ViQS0cknn5zYRQHV9J/vs+3bt9dDDz2ke+65R88++6xmz56tyy67TH/7299ILKqhQaIXgOo544wzVL9+fb3xxhuJXgoQaOzYserXr58uu+wyTZw4UQMHDtRBBx2U6GUB1bbv++zmzZt1++2366OPPtJ5552nrVu36sEHH9TMmTM1duxYEowQfHORJJo1a6bWrVuroKAg0UsBAmVnZ6tHjx6SpK5du+rYY4+t+DuQDPZ9n83Pz9fRRx+tN954Q8ccc4y6deumefPm6de//rU+//zzRC/1J49vLpIIt8cAQM0qf5/t37+/+vfvXyU+ZMiQ2l5SUiK5SBI7duzQli1b1KtXr0QvBQDqJNf77MiRIxOzoCTGr0WSxD//+U+VlpZygxwA1BDeZ/0huUgCGzZs0FVXXaUWLVro4osvTvRyAKDO4X3WL34t8hOzfPly7d27V3v37lVhYaHefPNNPfXUU6pfv77mzJmj9PT0RC8RAJIa77M1j+TiJ+aiiy6SJDVq1EgtW7ZU9+7d9ec//1mjRo3iggcAD3ifrXns0AkAALzingsAAOAVyQUAAPCK5AIAAHhFcgEAALwiuQAAAF6RXAAAAK8Sss9FWVmZCgoKlJqaqkgkkoglIMlFo1EVFxcrIyND9erVXo7MtYv9xbWLZFbd6zchyUVBQYEyMzMTMTXqmI0bN6p9+/a1Nh/XLnzh2kUyC7t+E5JcpKamSpKuvfZapaSkJGIJSHK7du3SHXfcUXEt1Zby+UaNGqVGjRrV6tyoG3bv3q3Jkycn7Np95JFH1KRJk1qdG3VHSUmJLrnkktDrNyHJRflXcikpKSQX2C+1/fVu+XyNGjVS48aNa3Vu1C2JunabNGmipk2b1urcqHvCrl9u6AQAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8IrkAgAAeEVyAQAAvCK5AAAAXpFcAAAAr0guAACAVyQXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVyQUAAPCK5AIAAHhFcgEAALwiuQAAAF6RXAAAAK9ILgAAgFcN4jlp0aJFWrRokQoLC1VWVlYpNmXKFC8LAwAAySnm5GL8+PG65ZZb1LdvX7Vr106RSKQm1gUAAJJUzMnFo48+qqlTp+q3v/1tTawHAAAkuZjvudi9e7eOP/74mlgLAACoA2JOLkaNGqUZM2bUxFoAAEAdUK1fi/zxj3+s+HNZWZkef/xxLVy4UIcffrgaNmxY6f+dNGmS3xXWks6dOyd6CT95+fn5iV4CAhx33HGJXsJP3pIlSxK9BARYsWJFopfwk9ejR49ELyEu1UouPvroo0p/7927tyRp+fLl3hcEAACSW7WSi7y8vJpeBwAAqCNivuciNzdXxcXFVY7v2LFDubm5XhYFAACSV8zJxbRp01RSUlLleElJiaZPn+5lUQAAIHlVe5+LoqIiRaNRRaNRFRcXKyUlpSJWWlqql19+WW3atKmRRQIAgORR7eSiZcuWikQiikQiysnJqRKPRCIaP36818UBAIDkU+3kIi8vT9FoVKeccopmzZqltLS0ilijRo2UlZWljIyMGlkkAABIHtVOLk466SRJ0tq1a9WhQwd6igAAgEAx9xZZv3691q9f74yfeOKJ+7UgAACQ3GJOLk4++eQqx/b9FqO0tHS/FgQAAJJbzKWoW7durfRfYWGh5s2bp6OPPlqvvvpqTawRAAAkkZi/uWjRokWVY4MHD1bjxo01ZswYLV261MvCAABAcor5mwuX9PR0rVq1ytdwAAAgScX8zcWyZcsq/T0ajWrz5s264447dMQRR3hbGAAASE4xJxe9e/dWJBJRNBqtdPzYY4/VlClTvC0MAAAkp5iTi7Vr11b6e7169ZSenl5pO3AAAHDgiumeiz179mjkyJH64YcflJWVpaysLGVmZpJYAACACjElFw0bNtTy5cvZnRMAADjFXC1y4YUX6sknn6yJtQAAgDog5nsudu/ercmTJ2vBggXq27evmjVrVik+adIkb4sDAADJJ+bkYvny5TryyCMlSatXr/a+IAAAkNxiTi7y8vJqYh0AAKCOiPmei9zcXBUXF1c5vmPHDuXm5npZFAAASF4xJxfTpk1TSUlJleMlJSWaPn26l0UBAIDkVe1fixQVFSkajSoajaq4uLjS3halpaV6+eWX1aZNmxpZJAAASB7VTi5atmypSCSiSCSinJycKvFIJKLx48d7XRwAAEg+1U4u8vLyFI1Gdcopp2jWrFlKS0uriDVq1EhZWVnKyMiokUUCAIDkUe3k4qSTTpL0Y2+RDh06sEsnAAAIFHMpalZWVk2sAwAA1BExV4sAAABYSC4AAIBXJBcAAMArkgsAAOBVtW7o7NOnT7WrQz788MP9WhAAAEhu1UoufvnLX1b8edeuXXr44YfVo0cPHXfccZKkd955R5999pkuvfTSGlkkAABIHtVKLm666aaKP48aNUpXXnmlJkyYUOX/2bhxo9/VAQCApBPzPRfPP/+8LrzwwirHL7jgAs2aNcvLogAAQPKKOblo0qSJFi9eXOX44sWLKzUzAwAAB6aYd+gcPXq0LrnkEi1dulTHHnuspB/vuZgyZYrGjRvnfYEAACC5xJxcXHvttercubPuv/9+zZgxQ5LUvXt3TZ06VcOHD/e+QAAAkFxiTi4kafjw4SQSAAAgUFzJhSTt3r1bhYWFKisrq3S8Q4cO+70oAACQvGJOLtasWaPc3Fy9/fbblY5Ho1FFIhGVlpZ6WxwAAEg+MScXI0eOVIMGDTR37ly1a9eu2jt3AgCAA0PMycXHH3+spUuXqlu3bjWxHgAAkORi3ueiR48e+vbbb2tiLQAAoA6IObm48847dc011+i1117Tli1bVFRUVOk/AABwYIv51yKDBg2SJJ166qmVjnNDJwAAkOJILvLy8mpiHQAAoI6IObk46aSTamIdAACgjog5uXjjjTfM+Iknnhj3YgAAQPKLObk4+eSTqxzbd68L7rkAAODAFnO1yNatWyv9V1hYqHnz5unoo4/Wq6++WhNrBAAASSTmby5atGhR5djgwYPVuHFjjRkzRkuXLvWyMAAAkJxi/ubCJT09XatWrfI1HAAASFIxf3OxbNmySn+PRqPavHmz7rjjDh1xxBHeFgYAAJJTzMlF7969FYlEFI1GKx0/9thjNWXKFG8LAwAAySnm5GLt2rWV/l6vXj2lp6crJSXF26IAAEDyijm5yMrKqol1AACAOiKuGzpff/11nXXWWerSpYu6du2qoUOH6s033/S9NgAAkIRiTi6efvppDRo0SE2bNtWVV16pyy+/XE2aNNGpp56qGTNm1MQaAQBAEon51yK33Xab7rrrLo0ZM6bi2P/8z/9o0qRJmjBhgn7zm994XSAAAEguMX9zkZ+fr7POOqvK8aFDh1a52RMAABx4Yk4uMjMztWjRoirHFy1apMzMTC+LAgAAySvmX4v86U9/0pVXXqmPP/5Yxx9/vCKRiBYvXqypU6fq/vvvr4k1AgCAJBJzcnHJJZfokEMO0T333KOZM2dKkrp3767nnntOw4YN875AAACQXGJKLvbu3avbbrtNubm5Wrx4cU2tCQAAJLGY7rlo0KCB7r77bpWWltbUegAAQJKL+YbOQYMG6bXXXquBpQAAgLog5nsuTj/9dF133XVavny5jjrqKDVr1qxSfOjQod4WBwAAkk9cN3RK0qRJk6rEIpEIvzIBAOAAF3NyUVZWVhPrAAAAdUTMyUVdlZ+fn+glAHFZsmRJopcAxKVHjx6JXgJqSLWTi5KSEi1atEhnnnmmJOm6667TDz/8UBGvX7++JkyYoJSUFP+rBAAASaPaycX06dM1d+7ciuTiwQcfVM+ePdWkSRNJ0sqVK5WRkVGpoRkAADjwVLsU9ZlnnlFubm6lYzNmzFBeXp7y8vJ09913V+zYCQAADlzVTi5Wr16tnJycir+npKSoXr3/P71fv35asWKF39UBAICkU+1fi3z//fdq0OD///dvvvmmUrysrKzSPRgAAODAVO1vLtq3b6/ly5c748uWLVP79u29LAoAACSvaicXZ5xxhsaNG6ddu3ZViZWUlGj8+PEaMmSI18UBAIDkU+1fi1x//fWaOXOmfvazn+nyyy9XTk6OIpGIVq5cqQcffFB79+7V9ddfX5NrBQAASaDayUXbtm319ttv65JLLtG1116raDQq6cctvwcPHqyHH35Ybdu2rbGFAgCA5BDTDp2dOnXSvHnz9N133+mLL76QJHXp0kVpaWk1sjgAAJB84tr+Oy0tTf369fO9FgAAUAdU+4ZOAACA6iC5AAAAXpFcAAAAr0guAACAVyQXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVyQUAAPCK5AIAAHhFcgEAALwiuQAAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8IrkAgAAeNUgEZNGo1FJ0q5duxIxPeqA8mun/FqqLeXz7d69u1bnRd1Rfu0k6totKSmp1XlRt5RfP2HXbyRa21e4pE2bNikzM7O2p0UdtHHjRrVv377W5uPahS9cu0hmYddvQpKLsrIyFRQUKDU1VZFIpLanRx0QjUZVXFysjIwM1atXe7/d49rF/uLaRTKr7vWbkOQCAADUXdzQCQAAvCK5AAAAXpFcAAAAr0guAACAVyQXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVyQUAAPCK5AIAAHhFcgEAALwiuQAAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8IrkAgAAeEVyAQAAvCK5AAAAXpFcAAAAr0guAACAVyQXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVyQUAAPCK5AIAAHhFcgEAALwiuQAAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8IrkAgAAeEVyAQAAvCK5AAAAXpFcAAAAr0guAACAVw0SMWlZWZkKCgqUmpqqSCSSiCUgyUWjURUXFysjI0P16tVejsy1i/3FtYtkVt3rNyHJRUFBgTIzMxMxNeqYjRs3qn379rU2H9cufOHaRTILu34TklykpqZKkkaPHq3GjRtXie/Zs8c8Pz093Rn76quvnLF27do5Y8XFxeacaWlpzlg0GnXGysrKzHFbtWrljG3ZssUZa9KkiTnuDz/84IzVr1/fGbPWG/YYtWjRwhnbuXOnMxb2r6iDDjqoyrFdu3bpuuuuq7iWakv5fFdddVXgtdupUyfz/C+++MIZ2759uzNmPUZB69jX8uXLnbGNGzc6Y6effro57gcffOCMZWVlOWNNmzY1x/3666+dsZYtWzpj69atc8a6du1qzrlixQpnLCcnxxkrKSkxx23UqFGVY7t379aMGTMSdu2+++67at68eZX4smXLzPPPOOMMZ+yGG25wxkaNGuWMTZo0yZzTeh97/PHHzXPjZX3+zJkzxzx3/vz5ztiTTz4Z13ruuOMOM26t13otDR482Bx32LBhgceLioqUmZkZev0mJLkof6Ns3Lhx4Btj2FeFKSkpzpj1RmudF5bQWB/m1gdyWHJhjWutNyy5sD6MGjRwP+2lpaXOWE09RmHJhTVubX+9u++1G/T8hH1wxnsN7k9y0bBhQ2fMSjTDxrWuo6AP1eqOa63XGtdaT039LNbrJezcRF27zZs3D/xgCLt2g5L8ctbja30IWY+PZP/DzVrP/rBeh2GPkXXtxrte6z1Dsj8vrcd3f55vKfz65YZOAADgFckFAADwiuQCAAB4RXIBAAC8SsgNneXKysoCb/QLu8nHqjzIzs52xrZu3eqMtWnTxpxz9+7dzli8N+JJ0oYNG5wx6ybIHTt2mOMG3Q1ezrqD2Lo5KOyGTuvO7rZt2zpjVpWEJBUWFlY5tmvXLvOcmvbtt9+G3iQYxLo5deHChc6YVbkxZcoUc86RI0c6Y126dHHGrGtIklnWeOyxxzpjmzZtMsf9/vvvnbFevXo5Y71793bGioqKzDmt94333nvPGbPeiyTp/PPPDzxn6tSp5nk1qX379oE36y1atMg875133nHG/vKXv8S1lieeeCKu8yTp3nvvdcbGjBljnjt58mRnzLopc8SIEea41nvCjTfe6IwNGTLEGbvsssvMOa2biq3PtYkTJ5rjzp49O/B42DVfjm8uAACAVyQXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABeJXSfi5SUlMCmLFYzJcmug4+3cUxY7a61P4G1v0Pr1q3Nca26aKvjotXQJ4zVidXaz+PQQw81x7XWtHnzZmcsbF+ToP019u7da55T01q0aBH43OXn55vnHX/88c5Y9+7dnbFt27Y5Y/379zfntPY1sfaymDVrljmutUeG1TE1IyPDHNfaE+XTTz91xqzXftieJIsXL3bGrD07whqXvf/++1WOWe8XteGbb74J3CcmbP+RsNe/yyuvvOKMrV271jzX2icobP8Hi9WpdcaMGc7Ys88+a45r7RN06623OmMXXHCBM2btySHZjc0WLFjgjPXp08cc19XYLKyRWjm+uQAAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8CqhpaglJSWBbcXDyjddJTLSj62wXawSGqvUVLJLjKxYWKvnoHbi5azSurS0NHNc63EIa6XtYpWpSnaJnVUKGFYGHFRCbD02taF58+aB19PWrVvN81atWuWM9evXzxlbuXKlM9asWTNzTqu023otHX300ea4AwcOdMas59TVyrmcVZJ78MEHO2NWmWpYebvV0j41NdUZW7dunTlu0DVilVfWhokTJwaWf5933nnmeRs2bHDG/vrXvzpjHTt2dMYuvfRSc86act999zlj1lYHVsmoJGVnZztjc+bMccasbRLCPkP+/ve/O2NWGfrIkSPNcV2vtbD1lOObCwAA4BXJBQAA8IrkAgAAeEVyAQAAvCK5AAAAXpFcAAAArxJaipqWlhZXV1RL586dnTGr619YKapVhmmVohYUFJjjNmjgfgqsNVlzSvF3DbVKCK1uqpK0fft2Z8zq6mmV+knB5Z2J7iy5atWqwHK+sM6RVgmt1Y3xhBNOcMa2bNlizmldY1999ZUzNmjQIHNcqzzWGjesjPiwww5zxoK6jJazujzOnz/fnLNbt25xxT7//HNz3KAOr0Hl97Vp+PDhgeXLU6dONc+zyk0nTJjgjHXt2rXaa4uF9di/8MIL5rlWl2arC671+SJJ//u//+uMWd2fR48e7Yy1adPGnNPqxHzFFVc4Y9b7giRNmjQp8HhQR90gfHMBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOBVQktRy8rKAsuyrK50kt0Zsbi4OK7zevbsac5pld9YJZh9+/Y1x12/fr0zZnUSDStna9mypTNmdWW0yp7efffduOe0OneGdfXcsWNHlWMlJSXmOTWtcePGgaVlVndDSfrss8+csd69eztjVnfD0047zZzTer6DSiXLhZVvWl1crU6YOTk55rhWmbHVRdgq9bNKeSW7rNYqnV2zZo05btB1GtZduKbNnTs38L0l7L3KKo2fOHGiM7Z48WJn7JprrjHn7NGjhzNmrff3v/+9Oa5VOvvll186Y9OmTTPHtco7X331VWesS5cuzljYdgbjxo0z4/H64x//GHi8qKhIN9xwQ+j5fHMBAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOBVtUtRW7VqZZa07eu7776Le0EAACC5VTu5uO+++yr+vGXLFt1666067bTTdNxxx0mSlixZovnz52vs2LHVnjw1NTWwdtqqV5ekaDTqjFmtxq39M6y9KiRpz549zpjVIt5qPS1JpaWlcY0b1Kp+X9ZeAU2bNnXGrPbd1r4bkr1ea1yr/bEkpaWlVTmW6L0CNmzYELj/QWZmpnmetY+IdS1Ye2AsW7bMnHPgwIHO2OrVq52x/v37m+O2bt3aGdu2bZszFnbtbtq0yRmzWtofeeSRzti8efPMOd98801nzGqzbe31IQW3iLfeo2rDmWeeGbi3zIABA+Ie02r7/Y9//MMZGzJkiDnukiVLnLFZs2Y5Yx9++KE57lFHHeWMnXvuuc7Y3LlzzXHbtWvnjFnX0VNPPeWM9erVy5wzNzfXjLu88sorZvz000+Pa9xy1U4uRowYUfHns88+W7fccosuv/zyimNXXnmlHnzwQS1cuFBjxozZr0UBAIDkFdc9F/Pnz9cvfvGLKsdPO+00LVy4cL8XBQAAkldcycXBBx+sOXPmVDn+j3/8w9xiGwAA1H1x9RYZP368fve73+m1116ruOfinXfe0bx58zR58mSvCwQAAMklruRi5MiR6t69ux544AHNnj1b0WhUPXr00FtvvaVjjjnG9xoBAEASibsr6jHHHKNnnnnG51oAAEAdEHdy8eWXX+qpp55Sfn6+7rvvPrVp00bz5s1TZmZmaPvyckVFRYElhWH7aVilXFYrcquEdefOneacVhmm1bbaKvWT7BbdVkmuVU4q2aWzVlmoNWdYibDVBt16/MJargeVoia65Xo0Gg281jZs2GCeZ/2sVpmqVdL885//3JzTugat9tKHH364OW69eu5btqwbu63XqGSXlB5yyCHOWFDZZ7nXX3/dnDPoBvVy1msp7B9YQa3eE11G3b59e6WmplY5Pnr0aPO8fbcj+E8LFixwxrKzs52xsD2R/vSnPzljVnlxWLlv8+bNnbGgEvNyYa3cZ86c6YyNGjXKGbvrrrucsXhLTSXphRdecMaGDRtmnvvII48EHq/ue29cN3S+/vrr6tWrl959913NmjWrYo+IZcuW6aabbopnSAAAUEfElVxce+21uvXWW7VgwYJK/5odOHCguekJAACo++JKLj799FP96le/qnI8PT3d/ModAADUfXElFy1btgzcsvmjjz4yt+cFAAB1X1zJxW9+8xv9+c9/1r///W9FIhGVlZXprbfe0lVXXaULL7zQ9xoBAEASiSu5uO2229ShQwcdeuih2r59u3r06KETTzxRxx9/vG688UbfawQAAEkkrlLUhg0b6plnntGECRP04YcfqqysTH369FHXrl1jGqdx48aBHRL3pyzU6ixpsTp6ho1rdSC1uuRJUmFhoTNmlcCFlbNZZbdW6axV8lhcXGzOaf2s1r04YaWJ69evr3Js165d5jk1rVmzZmbJmsvatWudMau82CrZKyoqMue0nm+rJO+LL76Ie9zhw4c7Y59//rk5bosWLZyxTz/91Bn7r//6L2csrDux1WHTKvUN+8dUUKn5zp079dxzz5nn1aR69eoFlhFbJcBhLr744rjOGzdunBl/+umnnTGrQeYpp5xijvviiy86Y0GtLcqFvdasbRisLqQNGrg/iletWmXOaXX0tcpfly9fbo7r2nqgup+xcX1zccstt2jnzp3q3LmzzjnnHA0fPlxdu3ZVSUmJbrnllniGBAAAdURcycX48eMr9rbY186dOzV+/Pj9XhQAAEhecSUX0Wg0cBfNTz75JHA3RQAAcOCI6Z6LVq1aKRKJKBKJKCcnp1KCUVpaqu3bt+sPf/iD90UCAIDkEVNycd999ykajSo3N1fjx4+vdONVo0aN1LFjx4oW7AAA4MAUU3IxYsQISVKnTp10/PHHx3W3PAAAqNviKkU96aSTKv5cUlJSpWTS6vQJAADqtriSi507d+qaa67RzJkzA/cvqG4d7N69ewP3ckhPTzfP+/rrr50xqw30pk2bnLGwPTqsb2m2bdvmjIXt4WDtZWG157bqosPiO3bsiGvOoJt497V161ZnzNoLJOwm4KB9TRLdct21V0B+fr55Xt++fZ2xoC31y5199tnOWNicWVlZzpjVmtp6ziTpm2++ccbefvttZ8zax0L6sY2AS79+/Zyx9957zxkLa8FtjWu9b4TtNxP0mgl7XGvatGnTAl9TYZV+N998szM2aNAgcz6XoH2O9mXtOWHtgWG1eQ+b1/rHsfV+LdnvZdYeGVYr97DfELz//vvOWEZGhjMW9ti79tap7h5DcVWLXH311frXv/6lhx9+WI0bN9bkyZM1fvx4ZWRkaPr06fEMCQAA6oi4vrl46aWXNH36dJ188snKzc3VCSecoC5duigrK0vPPPOMzj//fN/rBAAASSKuby6+++47derUSdKPXyGVf30yYMAAvfHGG/5WBwAAkk5cyUXnzp21bt06SVKPHj00c+ZMST9+o9GyZUtfawMAAEkoruTioosu0ieffCJJuu666yruvRgzZoyuvvpqrwsEAADJJa57LvbtRjdw4ECtXLlSH3zwgbKzs3XEEUd4WxwAAEg+cSUX/6lDhw7q0KFDzOfVr18/sNW51V5askt6yu8FCWK1crdKMCW7BKmgoMAZa9WqlTmuVQ5klfyElSdZ5XNdunRxxqx29mElj1aJnfXYW3NKwe3jE91yffv27YHlvt26dTPPs57voNby5d59911nLKwc0nq9nHPOOc7YihUrzHGt16lVLh5WenzNNdc4Y0888YQzZl1/Ye8pVnms9XrZuHGjOW7Y+0oiXH/99YHllmGvKevaHTBggDPWu3dvZ+zZZ58157zjjjucMavMctGiRea41vNilVF///335rhWyfONN97ojP33f/+3M5abm2vOOXLkSGds9uzZzlhY4YWrbX1Q09IgcScX7733nl577TUVFhZW2cth0qRJ8Q4LAACSXFzJxcSJE3XjjTfqZz/7mdq2bVtpo5iwjZYAAEDdFldycf/992vKlCnm1zEAAODAFFe1SL169dS/f3/fawEAAHVAXMnFmDFj9NBDD/leCwAAqAPi+rXIVVddpSFDhig7O1s9evSoUrlg3aEKAADqtriSiyuuuEJ5eXkaOHCgDj744Lhv4iwqKgoso2vatKl5XuvWrZ0xqwSuffv2zlhYiZFVptWkSRNnrF27dua4VhmrtdtpWPlmvN1NrZ/T6jgr2Z1Y27Rp44zF0zk2qDy1NuXk5AQ+B2FdE62SyKOOOsoZK9+0LoirZKyc9XqxyovD/pEwbNgwZ8wqTQ8rC/3nP//pjAWVrpezHvuwa8zqsvvVV185Y2vWrDHH7dixY5Vjie6KOmvWrMD32OLiYvM8q6S0fMfmIEGPQTmrk7IkffHFF86Y9ZxaXUYl6fbbb3fGhg8f7oydfvrp5rjvvPOOM/bvf//bGbNeo99++605p/U4nHXWWc6Y9XxK0tixYwOPV/f6jSu5mD59umbNmqUhQ4bEczoAAKjD4rrnIi0tTdnZ2b7XAgAA6oC4koubb75ZN910k7nrIgAAODDF9WuRBx54QF9++aXatm2rjh07Vrmh88MPP/SyOAAAkHziSi5++ctfel4GAACoK+JKLm666Sbf6wAAAHVEXPdcAAAAuFT7m4u0tDStXr1arVu3VqtWrcy9EsLq2MuVlJQE1uiG1aSHtWx2sdojd+/e3Tw3KyvLGbPqkMPa0x566KHOmLUHxo4dO8xxrZp/aw8Ma7+PsH0urL0nrP0JrFbOUnB7eWtfgtpwxBFHBO5vsnnzZvM8az8Aq6b/1FNPdcasa0iSVq1a5YwdfvjhzpjVjl2yr0/r9RS2p8wxxxzjjFk3kVvX9SOPPGLOae1H89FHHzljn332mTnuli1bqhyz2nLXhrPPPjuw5foVV1xhnmftV3HppZc6Y4cddpgzZu2Nsz/ntm3b1hz3hBNOcMbeeuutuGKS1KpVK2fsxRdfdMaszzRrvyNJWrhwoTN28803m+daJkyYEHi8qKhId999d+j51U4u7r33XqWmplb8me6nAAAgSLWTixEjRlT8mW6oAADAJa57LurXr6/CwsIqx7ds2WJ+/Q0AAOq+uJIL1+/Wf/jhBzVq1Gi/FgQAAJJbTKWoDzzwgKQfG19NnjxZzZs3r4iVlpbqjTfeULdu3fyuEAAAJJWYkot7771X0o/fXDz66KOVfgXSqFEjdezYUY8++qjfFQIAgKQSU3Kxdu1aSdLAgQM1e/Zss+ymOtLT0wPLEMNKE63ST6scbd9vWv5T2K9zrNbAVmldenq6Oa5VdmuV3fXs2dMc12qfvH79emfMaku/YcMGc854S1G//vprc9ygc63W8LUhPz8/sOX6+++/b56Xm5vrjP397393xqzHKKyMOqhktpx13Ye9JqzXkzVn0P1a+3rwwQedsVGjRjljVqt2qzRbstth9+3bN66YpMCqul27dpmtuWvaBx98EPjc/fWvfzXPe/bZZ50xq3zY2nQx7HGYO3euMxb0+is3Z84cc1yrbXifPn2csbPPPtscd+bMmc5YixYtnLGgkuVyRx55pDmn9d76t7/9zRmzHj/px3L7IGHbK5SL656LvLy8SolFaWmpPv74Y/MDGAAAHBjiSi5Gjx6tJ598UtKPicWJJ56oI488UpmZmXrttdd8rg8AACSZuJKL559/vuIrk5deeknr1q3TypUrNXr0aN1www1eFwgAAJJLXMnFli1bdMghh0iSXn75ZZ177rnKycnR7373O3366adeFwgAAJJLXMlF27ZttWLFCpWWlmrevHkaNGiQpB9vQGQTLQAADmxxtVy/6KKLNHz4cLVr106RSESDBw+WJL377rvscwEAwAEuruTi5ptv1mGHHaaNGzfq3HPPrShpqV+/vq677rpqj1NUVKTdu3dXOW6VCUl2iaZVXmM1WwsrbbTGtcruwn6WoJ+/XNeuXZ2xjz/+2BzX6mBqlZta3UbDShOtckmr6195QzyXoqKiKsdKS0vNc2ra559/HtitNejYvhYsWOCMWV17rU6anTt3Nuc86qijnDGrLNQqjZWkoUOHOmMrVqxwxoYNG2aOa5V35uXlOWPW6yVs/52xY8c6Y6tXr3bGwjqctmnTpsqxsO7CNe2VV14JfD8L62b99ttvO2NWGeuyZcucseeee86c0+raa71XrVy50hzX+kewtRXC9OnTzXE///xzZ+zkk092xqxS1DPPPNOcMycnxxm77777nLGwLR+uv/76wONhn2nlYrrKzzjjjIo9Hc455xzt3LmzUr30mWeeqWuvvTaWIQEAQB0TU3Ixf/78SlnLnXfeWSnb3bt3r1atWuVvdQAAIOnElFz85w6M1o6MAADgwJTYX/4BAIA6J6bkIhKJVLkp0rpJEgAAHHhiqhaJRqMaOXJkxZ3Gu3bt0h/+8IeKZmHVvYsUAADUXTElFyNGjKj09wsuuKDK/3PhhRdWe7zmzZsHlsOEbcQVVu7nYpWNhXX8tMpfO3Xq5IxZJXmSXaKZn5/vjIWVzlplRp988ol5rktYKWrLli2dMes5C/tZgr4dS/Q3ZikpKYGPh/UYSHbXXqskN6iksdy6devMOa1S1f98Te/r9ttvN8d94YUXnLFf/OIXzlhYx0qrdNH6WaxS3iFDhphzWqXdBQUFzljv3r3NcYM6SCa6o+/YsWN10EEHVTn+l7/8xTzv1ltvjWs+6/o899xzzXOt9+wTTzzRGbM660rSjh07nLFvvvnGGbPKksPGbdeunTNmldWGbTtglVn//ve/d8bOP/98c9yzzjor8PiOHTs0adIk81wpxuTiqaeeiuV/BwAAByBu6AQAAF6RXAAAAK9ILgAAgFckFwAAwCuSCwAA4BXJBQAA8Cquluu+uOq9w1oSWzXBrVq1csYaNHD/uFZMsvd4sNaTkZFhjmvV0Kenp5vnWqx9Oax9RKzH3qrhluyadOt5CdvX5OCDD65yLNF7BXTq1CmwbfUhhxxinmftg2E99tbz+e2335pzLly40Bn7+c9/7oxZ+01I0mGHHeaMffDBB85Y0B4L++rSpYszVlZW5oxZ6w3b4M/ar+KLL75wxsJeE0VFRTGvpaYtXrw4cL8Va58QSWrRokVc873zzjvO2MUXX2yea+3xYLU/79mzpznuggULnLGJEyc6Y8XFxea4t9xyizNm7e/y9ddfO2Nhj7u13ieeeMIZe+yxx8xxTzjhhMDjYe/X5fjmAgAAeEVyAQAAvCK5AAAAXpFcAAAAr0guAACAVyQXAADAq4SWou7ZsyewrKV58+bmeVbrWqucymp3HVbO17RpU2ds9+7d5rmW1q1bO2NbtmyJe06r3M8qu92zZ48zZj1+kl1iF41GnTGrzFIKbq+e6HK+b7/9NrA82Xr8JOn55593xoYPH+6MvfLKK85Y//79zTnffPNNZyyozLec9TqT7BJDq51zWNtqq7R727Ztzpj1vtGwYUNzzo8++sgZs0qs16xZY44b9PiWlpaa59S0pUuXKiUlpcpxq7RYssuLZ8+e7YxZpZIrVqww59y6daszZr135uXlmeNaz8HcuXOdMasUWpKGDh3qjFk/y/r1650xq4xfki644AJn7NJLL3XGzjvvPHPcDRs2BB4PK78uxzcXAADAK5ILAADgFckFAADwiuQCAAB4RXIBAAC8IrkAAABeJbQUNRqNBpb2hHUD/e6778wxXayuk1bXU8kuQbLKVMNKJq2yHmtN1pySlJqa6oxZj9/+/CxWyZT1s4SV1QY9p9bzXBuaNWsW2BX11VdfNc8LKgEsZ3Vc7NChgzNmPZ+S3YnV6h5plZpKUrdu3ZyxRx55xBkLK51dt26dM2ZdY/tz7VrX06ZNm5yx9u3bm+MGveeElSvXtOzs7MDHKuw1ZV2Dv/rVr5yxxx9/3BkrLCw057SetyFDhjhj//rXv8xxJ02aZMZdpkyZYsaDuuCWs0pGn3zySWdsf8ruH374YWfM6uAqSePGjQs8bv2M++KbCwAA4BXJBQAA8IrkAgAAeEVyAQAAvCK5AAAAXpFcAAAAr0guAACAVwnd5yISiQTWgVutaSW7ztZq12y19rb2hZCktLQ0ZyyobXw5q12zZNft79q1K67zJPsxsvbssPbdCGvBba3XGtdqAS8F13knuuV6dna2mjRpUuW41TpZknr16uWMffnll87YsGHDnLFp06aZc1qPvbUHweDBg81xX3zxRWfMWu9LL71kjnvKKac4Y3369HHG8vPz457TapXdqVMnZ6xjx47muEESfe0efPDBatasWZXjs2bNMs+z9rK48sornbFzzjnHGVuyZIk5Z9u2bZ2xxx57zBnLzc01x73rrrucsYsvvtgZs64/SXrrrbecMevxtfaqufrqq805rb03rMch7DPEtVeN9Rm7L765AAAAXpFcAAAAr0guAACAVyQXAADAK5ILAADgVUKqRcq777numraqDsLi8d6J3bBhQzNu3SFrVYuE/SxWJ0LrZwm7YzfexygSicQ9pzVuaWmpMxbWFTWoi2T5XLXdHbV8PtfjG9bx0nperMfBqnQKm9OKW5VDYc+L9ZzGey2EnWs9DtZja/2cYXNaj1887zflj2uirl1X9VDYz2JVoFmPkfWchV1j1pqsKrPt27eb41rXivVzho0bb7VivOsJmzPeykrJ/fiWrzXs+o1EE9C7etOmTcrMzKztaVEHbdy4MbTttU9cu/CFaxfJLOz6TUhyUVZWpoKCAqWmppr/UgZcotGoiouLlZGREbhXSk3h2sX+4tpFMqvu9ZuQ5AIAANRd3NAJAAC8IrkAAABekVwAAACvSC4AAIBXJBcAAMArkgsAAOAVyQUAAPDq/wCTuwzsANvf4AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,axes = plt.subplots(2,3)\n",
    "for par,plotpar,axs in zip(pars,['D*' if par == 'Dstar' else par for par in pars],axes.T):\n",
    "    estmap = ivim.io.base.read_im(outbase_fit+f'_{par}.nii.gz')\n",
    "    truemap = ivim.io.base.read_im(files_sim[par])\n",
    "    for im,ax in zip([truemap,estmap],axs):\n",
    "        ax.imshow(im[...,0], cmap='gray', vmin=0, vmax=1.5*np.max(truemap))\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        if ax == axs[0]:\n",
    "            ax.set_title(plotpar)\n",
    "    if axs[0] is axes[0][0]:\n",
    "        axs[0].set_ylabel('Ground truth')\n",
    "        axs[1].set_ylabel('Estimated')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ivim",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}