var_model_kimmel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 417,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import cv2\n",
    "\n",
    "from scipy.ndimage.filters import gaussian_filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_(s, k_pyr):\n",
    "    n, m = s.shape\n",
    "    s_new = np.zeros((n+2, m+2))\n",
    "    s_new[1:-1, 1:-1] = s\n",
    "    s_new[0, 1:-1] = s[0,:]\n",
    "    s_new[-1, 1:-1] = s[-1, :]\n",
    "    s_new[1:-1,0] = s[:,0]\n",
    "    s_new[1:-1,-1] = s[:,-1]\n",
    "    s_nxt = np.zeros((n, m))\n",
    "#     print(s_new)\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(m):\n",
    "            s_nxt[i,j] = np.sum(s_new[i:i+3,j:j+3] * k_pyr)\n",
    "#             print(i, j, s_nxt[i,j])\n",
    "            \n",
    "    return s_nxt\n",
    "\n",
    "\n",
    "def gauss(s, k_pyr):\n",
    "    n, m = s.shape\n",
    "    s_nxt = np.zeros((n-2, m-2))\n",
    "    \n",
    "    for i in range(n-2):\n",
    "        for j in range(m-2):\n",
    "            s_nxt[i,j] = np.sum(s[i:i+3,j:j+3] * k_pyr)\n",
    "            \n",
    "    return s_nxt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [],
   "source": [
    "def up_scale(l):\n",
    "    n, m = l.shape\n",
    "\n",
    "    ind_x = np.tile(np.arange(n), (m, 1)).T *2\n",
    "    ind_y = np.tile(np.arange(m), (n, 1)) *2\n",
    "    out = np.zeros((2*n, 2*m))\n",
    "\n",
    "    out[ind_x, ind_y] = l\n",
    "    out[ind_x+1, ind_y] = l\n",
    "    out[ind_x, ind_y+1] = l\n",
    "    out[ind_x+1, ind_y+1] = l\n",
    "\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decimate(s):\n",
    "    n, m  = s.shape\n",
    "    ind_x = np.tile(np.arange(0, n, 2), (int((m+1)/2),1)).T\n",
    "    ind_y = np.tile(np.arange(0, m, 2), (int((n+1)/2),1))\n",
    "#     print(ind_x)\n",
    "#     print(ind_y)\n",
    "    \n",
    "    return s[ind_x, ind_y]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_laplace(G, k):\n",
    "    return cv2.Laplacian(G, cv2.CV_64F) * 2**(-2*(k-1))\n",
    "\n",
    "def my_num_laplace(G):\n",
    "    kernel = np.array([[0,1,0], [1,-4,1], [0,1,0]])\n",
    "    out = gauss_(G, kernel)\n",
    "    return out\n",
    "    \n",
    "#     return cv2.Laplacian(G, ddepth=cv2.CV_64F,ksize=3) #* 2**(-2*(k-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {},
   "outputs": [],
   "source": [
    "# n, m = 100, 100\n",
    "# p = 3\n",
    "\n",
    "# k_pyr = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) /16\n",
    "# k = cv2.getGaussianKernel(3, 0)\n",
    "# k = k*k.T\n",
    "\n",
    "\n",
    "# print(k_pyr)\n",
    "# print(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [],
   "source": [
    "# s_nxt = gauss(s_ori, k_pyr)\n",
    "# s_nxt_ = gauss_(s_ori, k_pyr)\n",
    "\n",
    "# # s__ = cv2.filter2D(s, -1, k)\n",
    "# s__ = cv2.GaussianBlur(s_ori , (3, 3), 0)\n",
    "\n",
    "# print(s_nxt)\n",
    "# print()\n",
    "# print(s_nxt_)\n",
    "# print()\n",
    "# print(s__)\n",
    "\n",
    "\n",
    "\n",
    "# # s = s_nxt\n",
    "                 \n",
    "# print(decimate(s__))                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Main Loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [],
   "source": [
    "# k_lap = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]])\n",
    "\n",
    "# # Gb = gauss(s, k_lap)\n",
    "# Gb = cv2.Laplacian(s, cv2.CV_64F)\n",
    "\n",
    "# print(s0)\n",
    "# print(s1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x = .5\n",
    "# y = cv2.Laplacian(x, cv2.CV_64F)\n",
    "\n",
    "# print(y.shape)\n",
    "# print(x.shape, y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h2> Main"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 388,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = cv2.imread('C:/Users/User/Downloads/uni/CO421/image2.jpg', cv2.IMREAD_GRAYSCALE)\n",
    "cv2.imshow('im1', S)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 545,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0 255.99999999999991\n"
     ]
    }
   ],
   "source": [
    "# 0 to inf correction\n",
    "S_mod = S.astype(np.float64)\n",
    "S_mod += 1\n",
    "S_mod = np.where(S>255, 255, S_mod)\n",
    "\n",
    "print(np.min(S_mod), np.max(S_mod))\n",
    "\n",
    "cv2.imshow('im2', S_mod/255)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 546,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 (600, 1200)\n",
      "2 (300, 600)\n",
      "3 (150, 300)\n",
      "4 (75, 150)\n",
      "5 (38, 75)\n",
      "6 (19, 38)\n",
      "7 (10, 19)\n",
      "8 (5, 10)\n",
      "9 (3, 5)\n",
      "10 (2, 3)\n"
     ]
    }
   ],
   "source": [
    "s_ori = np.log(S_mod)\n",
    "\n",
    "# S_ = np.exp(s_ori)\n",
    "# cv2.imshow('im3', S_/255)\n",
    "# cv2.waitKey(0)\n",
    "# cv2.destroyAllWindows()\n",
    "\n",
    "\n",
    "p = 10\n",
    "\n",
    "s_pyr = {}\n",
    "\n",
    "s_pyr[1] = s_ori\n",
    "\n",
    "for i in range(1, p):\n",
    "    s = s_pyr[i]\n",
    "    s_nxt = cv2.GaussianBlur(s , (3, 3), 0)\n",
    "    s_pyr[i+1] = decimate(s_nxt)\n",
    "\n",
    "for i in s_pyr:\n",
    "    print(i, s_pyr[i].shape)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 547,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, p+1):\n",
    "    S_pyr = np.exp(s_pyr[i])\n",
    "    for i in range(i-1):\n",
    "        S_pyr = up_scale(S_pyr)\n",
    "    cv2.imshow('im3', S_pyr/255)\n",
    "    cv2.waitKey(0)\n",
    "    cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 548,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K 10 (2, 3) (2, 3) 4.191401326066419 4.191401326066419\n",
      "\n",
      "K 9 (3, 5) (4, 6) 4.295779620766383 4.1914547967736775\n",
      "\n",
      "K 8 (5, 10) (6, 10) 4.452170180092242 4.2959331780400545\n",
      "\n",
      "K 7 (10, 19) (10, 20) 4.878235675575677 4.452598776005438\n",
      "\n",
      "K 6 (19, 38) (20, 38) 5.2912480572332 4.890663311303346\n",
      "\n",
      "K 5 (38, 75) (38, 76) 5.534969970452513 5.328773532879069\n",
      "\n",
      "K 4 (75, 150) (76, 150) 5.541512941575129 5.588323001030166\n",
      "\n",
      "K 3 (150, 300) (150, 300) 5.542269835939125 5.667881925419334\n",
      "\n",
      "K 2 (300, 600) (300, 600) 5.544198969649278 5.749979968906675\n",
      "\n",
      "K 1 (600, 1200) (600, 1200) 5.545177444479562 5.801934655703318\n",
      "\n",
      "(600, 1200)\n"
     ]
    }
   ],
   "source": [
    "k = p\n",
    "\n",
    "a = .1\n",
    "b = 1000\n",
    "\n",
    "l = {}\n",
    "\n",
    "l[0] = s_pyr[p]\n",
    "# l[0] = np.array(np.max(s_pyr[p])).reshape((1, 1))\n",
    "\n",
    "Tk = 10\n",
    "\n",
    "\n",
    "while k > 0:\n",
    "\n",
    "    print(\"K\", k, s_pyr[k].shape, l[0].shape, np.max(s_pyr[k]), np.max(l[0]))\n",
    "\n",
    "    # print(s_pyr[k])\n",
    "    # UPSAMPLE error\n",
    "    if (s_pyr[k].shape[0] < l[0].shape[0]):\n",
    "        l[0] = l[0][:-1, :]\n",
    "    if (s_pyr[k].shape[1] < l[0].shape[1]):\n",
    "        l[0] = l[0][:, :-1]\n",
    "\n",
    "\n",
    "    Gb = num_laplace(s_pyr[k], k)\n",
    "\n",
    "#     print(\"K\", k, Gb.shape, s_pyr[k].shape, l[0].shape)\n",
    "    # print(l)\n",
    "\n",
    "    for j in range(1, Tk+1):\n",
    "        \n",
    "        Ga = num_laplace(l[j-1], k)\n",
    "        # print(Ga.shape, Gb.shape, s_pyr[k].shape)\n",
    "        \n",
    "        G = Ga + a*(l[j-1]-s_pyr[k]) - b*(Ga-Gb)\n",
    "        # print(G.shape)\n",
    "        # print(np.max(G),np.max(Ga), np.max(Gb))\n",
    "\n",
    "        u_a = np.sum(G*G)\n",
    "\n",
    "        u_b = -1* np.sum(G * num_laplace(G, k))\n",
    "\n",
    "        u_NSD = u_a/(a*u_a + (1+b)*u_b)\n",
    "#         print(\"u_NDS\", u_NSD)\n",
    "\n",
    "        l[j] = l[j-1] -u_NSD*G\n",
    "    #     print(l[j].shape)\n",
    "\n",
    "        l[j] = np.maximum(s_pyr[k], l[j])\n",
    "\n",
    "#         print (\"J\", j, Ga.shape, u_a.shape, u_b.shape, l[j].shape)\n",
    "    \n",
    "    out = l[Tk]\n",
    "    l = {}\n",
    "    l[0] = up_scale(out)\n",
    "    print()\n",
    "    \n",
    "    \n",
    "    k-= 1\n",
    "    \n",
    "print (out.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 549,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9204836461399304\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH7pJREFUeJzt3XmUZGWd5vHvL9bMjNwzI2vfAREpQExARQWlsQFpQW0dcBmw1TrS3W7TTg+2Z6ZnnOmZ7jNzuts5bk0jI7KqSCsiLogoshXUArVCUXtlbblvkUts7/wRkUWRVUVFZiw3IvL5nFMnI27civu7VVlPvfnGu5hzDhERqXw+rwsQEZHCUKCLiFQJBbqISJVQoIuIVAkFuohIlVCgi4hUCQW6iEiVUKCLiFQJBbqISJUIlPJi7e3tbvny5aW8pIhIxVu/fn2vcy56uvNKGujLly9n3bp1pbykiEjFM7N9uZynLhcRkSqhQBcRqRIKdBGRKqFAFxGpEgp0EZEqoUAXEakSCnQRkSqhQBcRqRIKdBGRKlHSmaIiIpXi3rX7jz3+6CVLPawkd2qhi4hUCQW6iEiVUKCLiFQJBbqISJVQoIuIVAkFuohIlVCgi4hUCQW6iEiVUKCLiFSJ0wa6md1hZt1mtuUkr33ZzJyZtRenPBERyVUuLfTvAVdNP2hmS4Argf3TXxMRkdI7baA7554A+k/y0j8Bfw24QhclIiIzN6vFuczs/cBB59yLZlbgkkREvHH8glyVaMaBbmZ1wFeB9+Z4/hpgDcDSpZWxYpmISCWazSiXVcAK4EUz2wssBjaY2fyTneycu8051+mc64xGo7OvVEREXteMW+jOuc1Ax9TzbKh3Oud6C1iXiIjMUC7DFu8DngHeYGZdZvap4pclIiIzddoWunPuxtO8vrxg1YiIyKxppqiISJVQoIuIVAkFuohIlVCgi4hUCQW6iEiVUKCLiFQJBbqISJVQoIuIVAkFuohIlVCgi4icxKHBcfb0xrwuY0YU6CIi02zYN8C3f7+Ln2w86HUpMzKrDS5ERKpR2jl+vfUoT7zSg89gLJ70uqQZUaCLiGT9cssRntzZy8UrWgkHfDy1sxfnKmeXTXW5iIhk7emNsbI9wnXnLyQSCpB2EE+lvS4rZwp0EZGs2GSSptogZkZt0A/AeDzlcVW5U6CLiGTF4kki4UxPdE0oG+gJBbqISEUZj6dIpByRbJAfa6Er0EVEKktfbBLgWAu9NhvsE+pyERGpLP2xOAB1oWygH2uhV9GHomZ2h5l1m9mW4479bzN7ycw2mdm/mVlzccsUESmuvmygR8LV3eXyPeCqacceBc51zp0H7AC+UuC6RERKauBYoGda6OGgD6PKRrk4554A+qcd+7VzbmoK1bPA4iLUJiJSMlNdLpFsl4vPjJqgv+pa6KfzZ8AvTvWima0xs3Vmtq6np6cAlxMRKby+WByfQU3w1VisDfmZmCuBbmZfBZLAPac6xzl3m3Ou0znXGY1G87mciEjRDMTiREIBzOzYsZqgr6K6XGa9louZ3QRcC1zhKmmxAxGRk+iLxanLfiA6pbbCulxmFehmdhXwn4DLnHNjhS1JRKT0+rMt9OPVBv0Mj096VNHM5TJs8T7gGeANZtZlZp8CvgE0AI+a2Qtm9p0i1ykiUlQDsfixES5TakNV1kJ3zt14ksPfLUItIiKe6YvF6WgMv+ZYbXAOfSgqIlINEqk0Q+OJY7NEp9QG/STTrmJCXYEuInPe4FgC4IQul6kVF4fGEyWvaTYU6CIy5706qejEUS6gQBcRqRjTV1qcokAXEakw06f9T5laQndoTIEuIlIRBqattDhFLXQRkQrTN20t9CkKdBGRCtMfi9NYE8Dvs9cc1ygXEZEK0x+L01YfPuG4z4xwwMfwhAJdRKQi9MfitEZCJ32tNuRXC11EpFL0x+K01J0i0IN+hhXoIiKVoT8Wp+1ULfSgWugiIhXBOcfAWJzW+pMHeo0CXUSkMgxPJEmkHK2n6nKpoD70We9YJCJSDaZmibZGQkwm0ye8fqoul3vX7j/2+KOXLC1egTOgFrqIzGnHAv0UXS6ZjaLTTCbLfwldBbqIzGlTgf56H4pCZUwuUqCLyJzWn11p8fWGLQIMjydLVtNs5bKn6B1m1m1mW4471mpmj5rZK9mvLcUtU0SkOPpjmZZ32+t0uUD1tNC/B1w17ditwGPOuTOBx7LPRUQqzlM7ewn6jZ9sPHTS119toVdBoDvnngD6px2+Drgz+/hO4PoC1yUiUhKxyeQJqywer2YO9KHPc84dBsh+7ShcSSIipTMWT52wDvrxqq3LJS9mtsbM1pnZup6enmJfTkRkRmLx5Ak7FR1vLoxyOWpmCwCyX7tPdaJz7jbnXKdzrjMajc7yciIixTE6maQ+fOpA9/uMugqZLTrbQH8IuCn7+Cbgp4UpR0SkdJxzxCaTJ2wOPV1TbbA6At3M7gOeAd5gZl1m9ing74ErzewV4MrscxGRijIWT5FIuddtoUMm0CthlMtp13Jxzt14ipeuKHAtIiIl1TuamVR0ukBvrJYWuohIteodzUz7nzNdLiIi1epYC73mNC30msroclGgi8ic1ZdtoZ+uy6WlLsjAmAJdRKRs9WVb6JHQqScWAbREQownUkwkynsJXQW6iMxZvaOT1AR9BPyvH4VTKzEOjMVLUdasKdBFZM7qjcVP290C0BoJAjAQK+9uFwW6iMxZvSOTOQV6c7aFPqgWuohIeeqLxU87ZBEy+40C9CvQRUTKU99obi30x1/KLFf1661Hi11SXhToIjInJVJpBsYSOQX61HrpY/Hy3oZOgS4ic9JALLdZopBZcTEc8DEW17BFEZGy05vjpKIpkXBAgS4iUo5yXZhrSl3Iry4XEZFy1BebeaDHJtVCFxEpO70j2S6X0yzMNaUuFFALXUSkHPXGJgn5fYQDucVgJORXH7qISDnqG43TXh/CzHI6vzYUYDKZJp5MF7my2VOgi8ic1Ds6SVt9OOfzI+HMiozlPP1fgS4ic1LfaJy2+lDO509NLirnddHzCnQz+5KZbTWzLWZ2n5nVFKowEZFi6hudpH0GLfS67Jrp/bEqbKGb2SLg80Cnc+5cwA/cUKjCRESKxTlH74xb6NXf5RIAas0sANQBh/IvSUSkuEYmk8RTaaIz6UPPdrmU84qLsw5059xB4P8A+4HDwJBz7teFKkxEpFh6RzKTimbXQq/CPnQzawGuA1YAC4GImX38JOetMbN1Zraup6dn9pWKiBRIX7YfvC2Sews94PcRCviqsw8d+CNgj3OuxzmXAB4E3j79JOfcbc65TudcZzQazeNyIiKFMbU59Ew+FIVMK72c9xXNJ9D3A281szrLjMy/AthemLJERIqnJ7vSYvsMulwg048+UI0tdOfcWuABYAOwOftetxWoLhGRoplqoU9tLZerTAu9fPvQc1uV5hScc38L/G2BahERKYm+0TgtdUEC/pm1aau5y0VEpCLNdNr/lLpq7XIREalUUwtzzVRd2M/wRJJkqjwX6FKgi8ick08LHWBwvDz70RXoIjKnpNKO/f1jDI0luHft/hn93kh2clG5drso0EVkTjk0OE4y7ehomH0LvVxHuijQRWRO2dkzCkB0VoFe3isuKtBFZE7Z1Z0N9Fn1oZf3iosKdBGZU3b1jFIX8lMXnvk0HHW5iIiUkV3dsVn1nwOEAj5qgr6ynVykQBeROWVXz+is+s+ntNaFNMpFRMRrA7E4fbH4rPrPpzTXhdRCFxHx2q48RrhMaY2E1IcuIuK1VwN99vvZN9cF1eUiIuK1XT0xwgEfzXXBWb9HpoWuQBcR8dSu7lFWtEfwmc36PToawgyMJZhIpApYWWEo0EVkztjZM8qqjvq83mP14mYADg6OF6KkglKgi8icMJFIcaB/jFXR/AL9gmygH+gfK0RZBaVAF5E5YV/fGGkHZ+TZQm+qC7KyPcKBgSproZtZs5k9YGYvmdl2M3tboQoTESmkqREuq6KRvN/rgiXNdPWP4ZzL+70KKd8W+teBXzrnzgbOB7bnX5KISOFNLcq1sj2/FjrABUubGZlMMlRmG13MOtDNrBF4F/BdAOdc3Dk3WKjCREQKaWfPKIuaa6nNrpiYjwuWZPvRy6zbJZ8W+kqgB/h/ZrbRzG43s/x/lhERKYJdPaN5959POXt+IwGfld0Ho/kEegC4EPi2c+7NQAy4dfpJZrbGzNaZ2bqenp48LiciMjujk0l2dhcu0EMBHwuba6sq0LuALufc2uzzB8gE/Gs4525zznU65zqj0WgelxMRmZ0HN3QxkUgT8vtmvI/oqSxpqeXg4DipdPl8MDrrQHfOHQEOmNkbsoeuALYVpCoRkQJJpx3fe3ovi1tqWdJaV7D3XdxaRzLtODI8UbD3zNfMt+x4rc8B95hZCNgNfDL/kkRECucPO3vZ3RPjI52LC/q+S1oy/zmUU7dLXoHunHsB6CxQLSIiBXfn03uJNoQ5d1FTQd+3pS5IJOSna6B8Al0zRUWkau3pjfHbl7r56MVLCfgKG3dmxpLWOg70l8/QRQW6iFSt7z+zl6Df+NglS4vy/otbaukZnWQsnizK+89Uvn3oIiJlZ19fjLuf3ce9a/fzvtUL6Gic/YYWr6e5NgRA70icpW3ex6n3FYiI5GEikWLDvgF29cbY2xtj++Fhnt7Vh99nXPWm+dx69RuLdu1IODPrtC82ydK2wo2gmS0FuohUnKGxBA+9eJDfvtTN07v6mEymAQj6jbZImPec3cHFy1v57OWrilpHJJyJ0P4y2ZJOgS4iFWNfX4xbf7yZ9fsGiKfStEVCXLishbM6GljQVENDTQA7bjeiQk0iOpVIKBOhfQp0EZHTS6Udv9/Rzd3P7ufxl7vxYZy/pIm3r2pnYXOtp7WphS4ikoNdPaP8dONBfrzhIAcHx4k2hPnLd59BJBygsWb2mzwXUijgI+g3BbqIyHRHhiZ4eNMhHnrxEJu6hjBgVbSeGy9eyjkLGvH7Zr+5c7FEQgH6RhXoIiIkUmke3nSI+587wHN7+3EOzl3UyDXnzue8xc001pZHa/xUIuEA/bFJr8sAFOgi4pHhiQT3rd3Pt363i6HxBO31Id5zdgfnL2qmvSHsdXk5i4T96nIRkbnp0OA4dzy5h7ue3cdkMs3KaITrL1jEmfPq8Vn5damcTiQUoGdULXQRmSOcc2zYP8CdT+/jkc2HccC5Cxt5x5lRFnk8UiVfkXCAl4+OeF0GoEAXkSL7mwc387uXuzk0NEFDOMBNb1/OJy9dzhM7er0urSAiIT9j8RQTiRQ1wfz3K82HAl1EiiKVdvzDL1/i3uf209EQ5roLFnLBkmbCAX9Fh/n0yUpTY9H7YnHPf9pQoItIwQ2NJfjc/Rt5YkcPl6xo5drzFpblkMNCODa5aFSBLiJVZl9fjA99+2kGYgk+cMEiLlrR6nVJRRUJZbpZestg6KICXUQKZsP+AT595zrG4yn+7B0rWNEe8bqkoju+he41BbqIFMQvNh/miz94gflNNdz89uW011fOWPJ8lNN6LnkHupn5gXXAQefctfmXJCKVZDye4u8e2cbdz+5nSUstH7tkGfXhudNWDGfXcymHFRcL8af+BWA70FiA9xKRCrKpa5Av/uAFdvfEeMcZ7Vx5zjyC/rm1s6WZ0RoJlcX0/7z+5M1sMfA+4PbClCMilcA5x+1/2M3133yKvtE4n3rHCq5ZvWDOhfmU1ki4Krpc/hn4a6ChALWISAUYHIvz5R9t4jfbj3LOgkY+dOFiakPeTqjxWlskVNldLmZ2LdDtnFtvZpe/znlrgDUAS5cWZ+dtESmNPb0xPvitpxgeT/K+1Qt4+6q21+wQNFe1RkIcGBjzuoy8ulwuBd5vZnuB+4H3mNnd009yzt3mnOt0znVGo9E8LiciXtrZPcq/+5dnmEymWfOulVx6RrvCPKs1EiqLYYuzDnTn3Fecc4udc8uBG4DfOuc+XrDKRKRs/OOjO7jum08xFk/xmXeuZEmr9zvcl5O2SIiRySSTyZSndczNTzBEJGc7u0e5/Q+78Rl8+p0rmNdY43VJZae1PgTAQCzhaR0FCXTn3O80Bl2k+hwdnuCmO57DzPjMO1fS0aAwP5m2SCbQ+zweuqgWuoic1PBEgpvueI7BsTg3v23uzPycjdZI5s/G66GLCnQROcFkMsX133iKHUdH+EjnEha1VPYmFMXWmm2hK9BFpKwkU2m+cN8L7O6N8cELF3PmPE0zOZ1jXS4ej3RRoIvIMem049YHN/PLrUd43+oFXLi0xeuSKkJTbRC/z9RCF5Hy4Jzjaw9v44H1XVxxdgeXntHudUkVw+czWuqCns8WnTtLoonIKaXSjv/2s618/5l9XLqqjfec3eF1SRWnHBboUqCLzHETiRRfuH8jv9p6lDXvWsmy1jrNAJ2FTKCry0VEPNIfi/Ox29fy661Hufa8BSxviyjMZ6ktElaXi4h4Y1PXILfcvYGe0UluuHgpqxc1eV1SRVMLXUQ88cPnD/Cn33kGgAc++zaFeQG0RkIMjiVIptKe1aAWusgckkil+drPtnHXs/tYFY1ww0VL2XJw2OuyqkLb1HouYwmiDd7MqlWgi8wRA7E4H/rO0+zuifHOM9p575vm4/epv7xQppZG6B2d9CzQ1eUiMgfs7B7h+m89xb6+Mf70wsVcvXqBwrzAOrIh3j3i3dBFtdBFqtzjL3fz+Xs3Eg76+cw7VrC0LeJ1SVVpqlXePTzhWQ0KdJEq5Zzju0/u4e9+vp35TTV84q3LaK4LeV1W1ZpaWrhnVC10ESmgiUSK//yTLfxofRdvWtjIh9+yhFBAPazFVBvy0xAO0D2sQBeRAjk0OM5n717Ppq4hPn/FmXQ0hPFpslBJRBvD9KgPXUQKYe3uPv78ng2MTib5+CXLmK/t4kqqoyFM94j60EUkT195cDM/XHeAlroQt1y+StvFeSDaUMPmrkHPrj/rTjUzW2Jmj5vZdjPbamZfKGRhIpK7u57Zy/3P7WdRcy2fvUx7f3ol00KvzC6XJPBXzrkNZtYArDezR51z2wpUm4ichnOOf/7NK3z9sVc4e34DN1y0VB9+eqijIcxYPMXoZJL6cOk7QGZ9RefcYeBw9vGImW0HFgEKdJEScM7xie8+x5M7e3nL0hauf/MiTRbyWEfjq2PR66P1Jb9+Qf4rN7PlwJuBtSd5bY2ZrTOzdT09PYW4nMicl047/stPt/Lkzl7eurKND1yoMC8H0fpMV5dX3S55B7qZ1QM/Br7onDthlR/n3G3OuU7nXGc0Gs33ciJz3lg8yX/44Qvc9ew+3nlmO39y3gINSywTUy10r4Yu5tXJY2ZBMmF+j3PuwcKUJCKn8vKREf7i3g3s6h7lynPmcflZUW1IUUa8Xs9l1oFume+i7wLbnXP/WLiSRGS6VNpx79p9/I+fb6ehJsjNly7nzI4Gr8uSaZpqg4T8Ps/GoufTQr8U+ASw2cxeyB77G+fcI/mXJSJTntjRw/98ZDsvHRnhjGg9H+5cTENN0Ouy5CTMjGhDmB6Ppv/nM8rlSUA/64kUgXOOtXv6+ebjO/nDK70saa3lhouWsHpRk7pYyly0IezZAl2aKSpSRhKpNL/ZdpR/eWI3LxwYpC7k5+pz5/O2lW0E/BpfXgk6GsLs6xvz5NoKdJEycKB/jB88f4AfrjtA98gkrZEQ7z9/IW9Z1kJQQV5ROhrDPL+335NrK9BFPBJPpnl021Huf34/T77SC8BZ8xr44zfN56x5DRpXXqGi9TUMjCWIJ9Mln7WrQBcpsQP9Y9z73H6+//ReYvEUTbVB3n12B53LWrQBRRWYGoveOzrJwubakl5bgS5SAslUmt/v6OHetfv57cvdGHD2/EYuXtHKGR31mhhURY4fi65AF6kiu3tG+beNB/nRui6ODE9QHw5w+VlRLlreqtZ4lZpa6dKLvUUV6CIFdqB/jEc2H+ahFw+x9dAwBpw5r54r3riUs+c3qm+8ykU9nC2qQBcpgAP9Yzy86TC/2HKYTV1DACxuqeWa1QtYvaiJplpNBJor2utDmCnQRSpK7+gkj2w+zE9fOMT6fQNAJsSvetN8zl3URGtEXSpzUcDvoy0S8mSBLgW6yAx0D0/w2EvdPLzpEM/s6iPtYH5jDX98zjxWL25WiAuQ2Yqux4P1XBToIqfgnONA/zhr9/Tx3J5+nt/bz97sDMC2SIjLzoqyelEz85u03Zu8VtSjregU6CLH6R6Z4KmdvfxhRy/P7O7j8FCmlVUX8rOsLcLV585nVbSeBU01WlNFTqmjIcyOIyMlv64CXeYk5xxHhifYcXSU7YeH2XZomK2HhtjVEwMyAb4qWs9Fy1tZ0R4h2hDWWHHJWUdDmN7RSdJph6+Eo5oU6FL1Eqk0Lx8ZYeuhIbYcHGbb4WF2HBlhZDJ57Jzm2iALmjJ94WfMa2BBU40CXGatoyFMMu0YGIvTVh8u2XUV6FJ1ukcmWL93gPX7Bth4YJAtB4eYTKYBCAd8LGiq4ZyFjcxrrKGjMcz8xhrqQvqnIIXT0Zj5XKVrYFyBLpKrobEELx0ZZvPBIV7sGuLFA4Ps7898cBnwGQuba+lc1sKS1joWNdfSEgmp5S1Fd9HyVkIBH/c/v5/zlzSX7LoKdKkIQ2MJdvaM8MrRUXYcHeWV7hF2HB3h6HE7wzTVBlnUXMvV585nWVuEhc01BHxaelZKL9oQ5sNvWcyP1nXxpT8661iLvdgU6OI55xzDE0mODk9wdHiCw4MTdA2Oc3BgnAP9Y+zqGaUvFj92ftCf2eZrYVMtb17SwrzGGhY212hbNikra961kvue288dT+3l1qvPLsk18wp0M7sK+DrgB253zv19QaqSquCcY3AswZFsUB8dnqB7eJLukUm6RyboGck87hmZPNbHPcWAhpoALZEQK9ojXLyilWh9mI7GGprrguo2kbK3rC3C1asXcM+z+/jzd6+isQQNjlkHupn5gW8CVwJdwPNm9pBzbluhipPSc86RSDniqTSJZJp4Kk08mSaRSpNMO+LJNJPJNJPJFJOJNMMTCUYmkgxPJOgbjWcDeoKjw5McHhpnIpE+4Rq1QT8NNQEaagK014dZ0RahviZAY22QxpogTbVBGmsD6i6RinfLZav4+abD3PPsfm65fFXRr5dPC/1iYKdzbjeAmd0PXAdUTaA7517z/PUmkjjnSDtIO4fLfk1POzZ1zrGvTB3PPE47SKdd9v2yv4/XvmcqnfmVdo5k2pFOZ74mU45kOk0q+zyVdiRSaSYTmfAdi6cYnUweC9/h8QRD4wmGx5PE4klik0nG4iniqTTTbjtnIb+PhpoA9eFMOHcua82Gc5DGmgANNUEaagLaUk3mjHMXNfHOM9u546k9fPLS5dQE/UW9Xj6Bvgg4cNzzLuCS/Mo5uf/+8Dbuf27/sefH581UGL76+LUnTYXm9N+XOf/VI+7YsUJVXX4MCAd91AT91GZ/1QT9ROvDLGquJeT3EfAbfp+PgM+yj42Az4ffB36fD79lFh8K+I2gz0c44KMm5Kcm4C/5dlsileCWy1bx0dvX8pvtR7n2vIVFvVY+gX6y5uoJcWhma4A12aejZvZyju/fDvTOsrZKpPutbrrfKvaxHO73T/4hr0ssy+WkfAK9C1hy3PPFwKHpJznnbgNum+mbm9k651zn7MurLLrf6qb7rW7lcr/5/Iz8PHCmma0wsxBwA/BQYcoSEZGZmnUL3TmXNLO/BH5FZtjiHc65rQWrTEREZiSvcejOuUeARwpUy3Qz7qapcLrf6qb7rW5lcb82fWieiIhUJo0zExGpEp4HupldZWYvm9lOM7v1JK/fbGY9ZvZC9tenvaizUE53v9lzPmJm28xsq5ndW+oaCymHv99/Ou7vdoeZDXpRZ6HkcL9LzexxM9toZpvM7Bov6iyUHO53mZk9lr3X35nZYi/qLBQzu8PMus1syyleNzP7v9k/j01mdmFJC3TOefaLzIepu4CVQAh4EThn2jk3A9/wss4S3++ZwEagJfu8w+u6i3m/087/HJkP1z2vvYh/v7cBt2QfnwPs9bruIt/vj4Cbso/fA9zldd153vO7gAuBLad4/RrgF2Tm6bwVWFvK+rxuoR9bPsA5Fwemlg+oVrnc72eAbzrnBgCcc90lrrGQZvr3eyNwX0kqK45c7tcBjdnHTZxk7kYFyeV+zwEeyz5+/CSvVxTn3BNA/+ucch3wfZfxLNBsZgtKU533XS4nWz5g0UnO+1D2x5cHzGzJSV6vFLnc71nAWWb2lJk9m13RslLl+veLmS0DVgC/LUFdxZLL/f5X4ONm1kVmhNjnSlNaUeRyvy8CH8o+/gDQYGZtJajNKzl/zxeD14Gey/IBPwOWO+fOA34D3Fn0qoonl/sNkOl2uZxMi/V2MyvdlieFldPyEFk3AA8451JFrKfYcrnfG4HvOecWk/nx/C4z8/rf4Wzlcr9fBi4zs43AZcBBIHnC76oeM/meLzivv5FOu3yAc67POTe1Lc2/Am8pUW3FkMtyCV3AT51zCefcHuBlMgFfiXJaHiLrBiq7uwVyu99PAT8EcM49A9SQWQekEuXy7/eQc+6Dzrk3A1/NHhsqXYklN5Pv+YLzOtBPu3zAtP6n9wPbS1hfoeWyXMJPgHcDmFk7mS6Y3SWtsnByWh7CzN4AtADPlLi+QsvlfvcDVwCY2RvJBHpPSassnFz+/bYf9xPIV4A7SlxjqT0E/PvsaJe3AkPOucOlurinW9C5UywfYGZfA9Y55x4CPm9m7yfzY1o/mVEvFSnH+/0V8F4z2wakgP/onOvzrurZy/F+IdMNcb/LDhOoVDne718B/2pmXyLzo/jNlXrfOd7v5cD/MjMHPAH8hWcFF4CZ3Ufmntqzn4P8LRAEcM59h8znItcAO4Ex4JMlra9Cv5dERGQar7tcRESkQBToIiJVQoEuIlIlFOgiIlVCgS4iUiUU6CIiVUKBLiJSJRToIiJV4v8DzLIIR3rfE1kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ill = out\n",
    "\n",
    "ref = s_ori - ill\n",
    "\n",
    "R = np.exp(ref)\n",
    "\n",
    "sns.distplot(R.flatten(),  bins=100)\n",
    "print(np.mean(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 550,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv2.imshow('ref', R)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 551,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.distplot(ill.flatten(),  bins=100, hist_kws=dict(cumulative=True), kde_kws=dict(cumulative=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 552,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2a0256ead30>"
      ]
     },
     "execution_count": 552,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmU3GWd7/H3t6qrl+p97056yx7CDjEswYAIioyCzjgjMO7O4DgyV69eHb13jjp6z9WjzuiM24i4IwiCIqMgKouEJSuBBEhn6046naTXJJ3el6rn/tHVsQmdpDpdVb9aPq9zcuiq/uVX3zp0Pv3U8/v+nsecc4iISHrxeV2AiIjEnsJdRCQNKdxFRNKQwl1EJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQNKdxFRNJQllcvXFFR4Zqamrx6eRGRlLR58+Ye51zl6Y7zLNybmprYtGmTVy8vIpKSzGxfNMdpWkZEJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQNKdxFRNKQwl1EJA0p3EVE0pDCXUQkDXl2h6rIXNy1vu3417dc0uBhJSLJSSN3EZE0pHAXEUlDCncRkTSkOXdJeZp/F3k1jdxFRNKQwl1EJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQNKdxFRNKQwl1EJA0p3EVE0lBU4W5m15nZDjPbbWafmuH7DWb2uJltMbOtZnZ97EsVEZFonTbczcwPfAt4E7ACuNnMVpxw2L8A9zrnLgRuAr4d60JFRCR60YzcVwG7nXMtzrkx4OfAjScc44CiyNfFwMHYlSgiIrMVTbjPB/ZPe9weeW66zwHvNLN24CHgn2Y6kZndamabzGxTd3f3GZQrIiLRiCbcbYbn3AmPbwZ+5JyrA64Hfmpmrzq3c+5259xK59zKysrK2VcrIiJRiSbc24H6aY/rePW0yweAewGcc88CuUBFLAoUEZHZiybcNwJLzGyBmWUzecH0wROOaQNeD2BmZzEZ7pp3ERHxyGnD3Tk3AdwGPAJsZ7Ir5iUz+7yZ3RA57OPA35vZC8DdwHudcydO3YiISIJEtROTc+4hJi+UTn/uM9O+fhlYHdvSRETkTOkOVUk5R4fGeHR7J13HRrwuRSRpaQ9VSTl3b9jPo81dPNbcxYp5RVy5tJK60qDXZYkkFY3cJaXctb6Nezbup6owh6uWVbKne4DvPLGHQ33DXpcmklQU7pJShsdCtB0eZEVtEdeuqOEjr1+KA3Z3DXhdmkhSUbhLStnV1U/YwbKaQgCK8wKUBgPsPzzkcWUiyUXhLillR0c/eQE/9WV/nmOvLwuy/4imZUSmU7hLygiHHTs7+1laXYDP/rwqRkNZkL7hcfqGxz2sTiS5KNwlZWw70MfgWOj4lMyU+kinTJumZkSOU7hLynisuQsDllS9MtxrS3LJ8pnm3UWmUbhLynhiRxf1ZUHyc155e0aWz8e8kjyFu8g0CndJCd39o7zQ3sfS6sIZv99QFuTA0WHGJsIJrkwkOSncJSU8tXtykdET59un1JcFmQg7th86lsiyRJKWwl1SwssHj5GT5aO2OHfG7zdEWiO3tB1JZFkiSUvhLilhV9cAiypf2QI5XXFegKLcLLbsP5rgykSSk8JdUsKuzgGWVBec8pj6siBb2hTuIqBwlxQwODrBgaPDLKk6dbg3lAVpOzxEz8BogioTSV4Kd0l6U4uCLa6a+WLqlKmbmba2a/QuonCXpLcrEu5LTzMtU1WYA0BL92DcaxJJdgp3SXq7uvrJ9vuOd8ScTF62n6LcLPb16mYmEYW7JL3dnQMsrMwny3/qH1czo7E8n326U1VE4S7Jb1fXAItPczF1SmN5kH29mpYRUbhLUhseC7H/yNBJlx04UWN5kANHhhkPaRkCyWwKd0lqe7oHcA46+ka4a33baY9vLM9nIuw4eFSbd0hmU7hLUtvV1Q/8uRPmdJrK8wF0UVUynsJdktquzgF8BuUF0YV7Y/lkR43m3SXTKdwlqe3sHKCiIAe/b+Y1ZU5UVZhDbsCnkbtkPIW7JLXdXf1UFc28EuRMzIzGsnz2KtwlwyncJWmNjIdoOzwU9Xz7lMbyIG2HNS0jmU3hLkmrpXuQsIv+YuqUyV73IcJhF6fKRJKfwl2S1vFOmVlMy9y1vo2u/lFGJ8J09Wt1SMlcCndJWnu6JjtlKgqyZ/X3yvInj9+rjhnJYAp3SVqtvUPUlQbJ8s3ux7Q8f3Iap00XVSWDKdwlae3tGaSpIn/Wf684L4DPNHKXzKZwl6TknGNv7yBN5ade5ncmfp9RGszW6pCS0RTukpQOD47RPzJxfDmB2SovyNZdqpLRFO6SlKZuQmqqmP3IHaAsP4d9vUM4p3ZIyUwKd0lKe3smR91nPHLPz6Z/ZIIjQ+OxLEskZSjcJSnt7R3EZ1BXemYj9/JIO6SmZiRTRRXuZnadme0ws91m9qmTHPM3Zvaymb1kZnfFtkzJNHsjbZDZWWc2/ig7Hu66qCqZKet0B5iZH/gWcC3QDmw0swedcy9PO2YJ8GlgtXPuiJlVxatgyQxn2gY5pTQ/GzOFu2SuaIZFq4DdzrkW59wY8HPgxhOO+XvgW865IwDOua7YlimZZC5tkFMCfh81Rbns0wJikqGiCff5wP5pj9sjz023FFhqZk+b2Tozuy5WBUrmmWsb5JSGsqBG7pKxogn3mXZJOLG/LAtYAlwF3AzcYWYlrzqR2a1mtsnMNnV3d8+2VskQU3eW7u0djGrf1JOZWh1SJBNFE+7tQP20x3XAwRmO+bVzbtw51wrsYDLsX8E5d7tzbqVzbmVlZeWZ1ixpbm/PZCBX5M9uqd8TNZbn0zMwyuDoRCzKEkkp0YT7RmCJmS0ws2zgJuDBE455AHgdgJlVMDlN0xLLQiVz7O0dxICS/MCcztNQNjln36ZlCCQDnTbcnXMTwG3AI8B24F7n3Etm9nkzuyFy2CNAr5m9DDwOfMI51xuvoiW97e0dojQ/e9arQZ7oz5tlK9wl85y2FRLAOfcQ8NAJz31m2tcO+Fjkj8ic7O0ZPH4T0lw0lk1ekNWWe5KJdIeqJJWpNsjyWW7QMZPiYICSYEAjd8lICndJKlNtkOVzvJg6pbEsqDl3yUhRTcuIJMpUG2QsRu4ADeX5vLD/aEzOlSjT2z9vuaTBw0oklSncJanEqg1ySmNZkIe2HWI8FCbgT54Pqu1HhvjGo7sJ+H0U5mbx7ssbqSnKxWym20pEZk/hLkllajXIubZBTmkoDxIKOw4cGZ7TWjWxMDIe4gM/2khzRz+9g2Ov+N53/rSHVU1l/POblnlUnaQbhbskldaewTPaFPtkGiO97vsOD3ke7j9+Zi9P7+llWXUhly4sZ0FFPs5B/8g480vz+N7aVv7qO89yVk0hbzl/HiXB2ExNSWZSuEtS2ds7yIIYhnBjZH2att5BwLu7okfGQ3xvbSuLKwt4z+VNJ3w3j1suaeBdlzXyw6f38h+P7uLuDW188MpFXpQqaSJ5JiEl4znnaO2ObbhXFeaQk+XzvB3yno376RkY5arlJ/8FE8zO4sOvW8zbLpjP/iPDPLWrJ4EVSrrRyF2SRvfAKINjoZiGu89nk6tDetgOOTYR5r/+tIeVjaUsOMlKl9M7ZM6rK2bbgT7+uL2T3V0DLK4qSFSpkkY0cpek0do92QYZy3CHyWUI2jwcuf9qSzuH+kb48NWLo+qGMTNuvGAeAb+PT9z3AqGwNvmW2dPIXZLGVI/7gop82o8Mz/l8U6Ph4bEQbYeHcM4ltNXwrvVthMKOr/1xJ/NKcjl4ZDjq1y/MDXDD+fO4Z9N+fvh0K3/32oVxrlbSjUbukjRaegbJ9vuYV5IX0/OW5WczPB6iu380pueNxq7Ofg4PjnHV0qpZ/2I5r66YJVUF/Nvvd/L9ta1xqlDSlcJdksbenkEayoP4fbEdXZdFbojyYt59R2c/2Vk+ltcWzvrvmhlvPLuG4fEQa3dpcxuZHYW7JI3Wnth2ykyZWmEy0R0zzjl2dvazqLLgjPv255XkcV5dMU/v6aGrfyTGFUo6U7hLUgiHHXt7h+IS7iX5AbJ8xp7ugZif+1R6BsY4MjTO0uq5dbtce1Y1obDjm4/tjlFlkgkU7pIUDvYNMzYRjku4Z/l8LK4qYEdHf8zPfSo7Oydfb0nV7KdkpisvyGFlYxl3rW/ztOtHUovCXZLC1IJhTSfpA5+r5TWFNB86Fpdzn8zOzn4qCnIoi8HGI1cvrwLgo/dsmdOm4ZI5FO6SFFp7JqdMXth/NC7htby2iIN9I/QNjcf83DMZGQ/R2jPIsjlOyUwpyguwsqmMF/b3cWwkMe9BUpvCXZJCa88Q2ZHlb+Nhec3k1EhzR2JG7+taepkIO5ZUz21KZrrVi8oJO8f6Fm1PLKencJek0NozQHlBdtxuMmo+NDn/fee6fXE5/4me2NFNls9ieg2hvCCH5TWFrG89zMh4KGbnlfSkcJeksLd3iPKC2GzQMZPC3CyC2X46jiWmnfDJnd0srMyP+QYhly+uYGgsxANbDsT0vJJ+FO7iufFQmLbDQ1TEaGu9mZgZNUW5dPTFP9zbeodo6RlkaQynZKYsrMintjiXHzzdinNac0ZOTuEunms/Mkwo7GK2td7J1BTn0nFshHCcF+J6rLkTgKVzbIGciZlx+aIKdnYO8NRuLQksJ6dwF89NdcrEc+QOUFOUy3jI0RbnZQge2tbBsupCKgrj88vq/LpiKgpyuP3JlricX9KDwl081zq1KXYc59xhcuQO8e2Y6Tw2wsZ9h7n+3Nq4vUaW38f7r2hi7a4ent9/NG6vI6lN4S6ea+0ZoDgvQDAnvitQVxXmYsD2Q/G7U/XhbYdwDv7ivJq4vQbAuy9rojgvwDcf2xXX15HUpXAXz7V0D7KwMv6bV2dn+SgvyInryH1qSmZxHObbpyvIyeL9qxfwx+1dvHSwL66vJalJ4S6e29M9wMKKxGwlV1OcS3Oc1pjpSsCUzHTvXd1EYU6WFhSTGWknJvHUwOgEncdGWVQV/5E7TF5UffFAH4OjE+THeBro4Rc7cG5yqd94r/8ydf6Lm0p5+MUOvvaHnVQX5XLLJQ1xfV1JHRq5i6daIsvwJmrkXhu5qLqjM/aj999uPURVYQ5VRbkxP/fJrF5UQbbfx+M7uhL2mpIaFO7iqZbIptiLEzhyhz8vRxArU1My584vjul5Tyc/J4tLFpaxrb2P3oHEbyMoyUvhLp7a0z2A32c0lCUm3EuCAYrzAmxpOxLT8/76+YM4B+ckONwBVi+uwO8zntRWfDKNwl081dI9SENZkOysxPwomhlrllby+I7umN2pOjIe4va1LVy6sIzqBE7JTCnKDXBxYynP7TuakOUVJDUo3MVTe7oHWJSANsjp8gI+egZG+ervd8TkfD9b30Z3/ygfvWZpTM53Jl67pBKH43trddeqTFK4i2dCYUdLzyALKxNzMXXK0qpCDGLSEjkyHuLrf9jJwor849cPvFCWn835dSXctb6Nw4NjntUhyUPhLp45eHRy39REj9yDOVk0lAVjcjPTnev20T86wevPqo5BZXOzZmklw+MhPn7vC16XIklA4S6e2R1pg9zTNZjwfUGX1xRy8OgInXNY3314LMR//amFhZX5cdnYe7aqi3JZUVvEupZeRie0mUemU7iLZ6amMSrjtHriqSyrLQLg8eYz7w//wdOt9AyM8vrl3o/ap6xaUMbweGhO70vSQ1ThbmbXmdkOM9ttZp86xXFvNzNnZitjV6Kkqz3dA+QF/DG/UzQa1YU5lAQDPHqGIfjVR3bw77/fydnzipJi1D5lcVUBhblZ3Le53etSxGOnDXcz8wPfAt4ErABuNrMVMxxXCPwPYH2si5T0tKdrwJNRO0y2RC6rLuSpXT2z3o90eCzEPZv2E8zx87YL5sepwjPjM+OC+hKe2NFNj25qymjRjNxXAbudcy3OuTHg58CNMxz3BeDLgBptJSotPYNUxnkN91NZXlPE8HiIdS29s/p7X3x4O939o7z94rq4L1N8Ji5qKGUi7Pj18we9LkU8FE24zwf2T3vcHnnuODO7EKh3zv0mhrVJGjs2Mk53/6hnI3eAhZX55AX8s9ps+tHtnfzk2X2sXlTOkjgv63umqotyOa+uWFMzGS6acLcZnjt+a5+Z+YCvAR8/7YnMbjWzTWa2qbtbt0pnMi8vpk4J+H285/ImHnj+IM9EsR/p5n1HuO2uLZw9r4g3nB3fzTjm6u0X17H90DGt9Z7Bogn3dqB+2uM6YPrnvULgHOAJM9sLXAo8ONNFVefc7c65lc65lZWVlWdetaS8PV1T+6Z6F+4AH71mCU3lQT79q20Mj5187v3lg8d43w83UF2Uw4/et4qAP7kbzd5y3jwCfuP+zdF/KpH0Es1P6EZgiZktMLNs4CbgwalvOuf6nHMVzrkm51wTsA64wTm3KS4VS1rY0z1Als8oy4/vptinkxvw8//+8lz29Q7x9Ud3znjMnu4B3v2D9eTnZHHn313i6aeNaD38YgdLqwu5Z9N+QjFaQ0dSy2nD3Tk3AdwGPAJsB+51zr1kZp83sxviXaCkp5buQRrKg/h9M836Jdbliyq46TX13LG29RWrRTrn+Nn6fdzwjadwDm56TQNP7uxJ+A1XZ+rsecUMjk6wtV2baGeiqC71O+ceAh464bnPnOTYq+ZelqS7XV39LKlK7Joyp/Lp68/it9sO8bZvP8P59SW8YUU1T+/u4Zk9vSyqzOcvL6yj1ONPGbO1tKoAA57Y0c2FDaVelyMJltwTh5KWxibC7O0dSqpuk+K8AB+6chFvWFGNc46vPLKDre19vPWC+bx/9YKUC3aYXEOnvizIE9qlKSMlX5OupL29vYOEwo4l1QUMjnq7Bsr0KZaSYDZXLavilksa6OgbITfg46FtHR5WN3dLqwv54/ZOz9tOJfE0cpeE29U52SmzOImmZaa7a30bjzV3pXywAyyrmfx09OROtR5nGoW7JNyurn7MYFGC13HPRLXFuRTkZPHjZ/emzIVgiQ2FuyTco9u7KA1m88vn1IMdbz4zllYXsqtzQC2RGUbhLgnX1T9CleZ/E2ZZTSHD4yHajwx5XYokkMJdEmoiFKanf4yqwsRvJJ2pFlcW4DPYEYNtBSV1KNwlofYdHiLkHFVFGrknSl62n4ayfHZ2KtwzicJdEmqqU0bTMom1tLqAg30jWuM9gyjcJaF2d02OHtVznVgLI7tFbWw97HElkigKd0moXV0DlAQD5GT5vS4lo8wrzSPgN9Yr3DOGwl0SalfngKZkPJDl89FQFlS4ZxCFuyRMKOzY0z2gThmPNFXk09xxjL6hca9LkQRQuEvC7D88xOhEWCN3jyyoyMc52LRPo/dMoHCXhNkV2X2pqkgjdy/UlwbJ9vs0NZMhFO6SMLsinTIauXsj4PdxQX2Jwj1DKNwlYXZ3DlBTlEtuQJ0yXlm1oIwXD/QxMDrhdSkSZwp3SZidXf0sqdZKkF5ataCMUNjx3L4jpz9YUprCXRLip8/uo/lQP04LE3rq4sZS/D5jfWuv16VInCncJSF6B0eZCDtqdDHVU/k5WZwzv5gNmndPewp3SYjOY5NrmtQUK9y9dNf6Nopys3iu7Sgj495ucSjxpXCXhOjoG8ZnWlMmGSyqLCAUdjzboqmZdKZwl4ToODZKeUEOAb9+5Ly2oCKfgN94vLnL61IkjvQvTRKio29Y8+1JIuD3saiygEe3d+F0hTttKdwl7gZGJzgyNK759iSyvKaIA0eHj981LOlH4S5xN7W9m0buyWNZTSEAj2lqJm0p3CXumjuOAQr3ZFKcF+Cs2iKFexpTuEvc7ejoJyfLR0kw4HUpMs3VyyvZvO+IlgBOUwp3ibvmQ/3UFOViZl6XItNcvbyKUNjxp13dXpcicaBwl7hyzrG94xjVupiadC6oL6U0GFBLZJpSuEtcHeoboX9kQvPtSeiejftpLM/nkZc6CIXVEpluFO4SV7qYmtyW1xQyNBZive5WTTsKd4mr5qk2SE3LJKWzaovIC/i5e+N+r0uRGFO4S1w1H+pnfkmeNuhIUgG/jwsbSvjdi4foHRj1uhyJIYW7xM1d69tY19JLYW6W16XIKbymqYzxkOP+59q9LkViSOEucTM2Eaa7f5Ta4jyvS5FTqC7K5TVNpdy9Yb/WmkkjCneJm46+YRwwv0ThnuxuXtVAa8+glgFOIwp3iZsDR4cBmF+qcE92/SMT5AX8fOnhZq9LkRhRuEvcHDg6Qn62nyLNuSe9qQurLx04Ro8urKaFqMLdzK4zsx1mttvMPjXD9z9mZi+b2VYze9TMGmNfqqSag0eHmV+ap2UHUsQlC8oJO8d3/7TH61IkBk4b7mbmB74FvAlYAdxsZitOOGwLsNI5dx5wH/DlWBcqqWVkPERX/wjzNN+eMioLc7iooZQfP7Pv+JSapK5oRu6rgN3OuRbn3Bjwc+DG6Qc45x53zg1FHq4D6mJbpqSa7YeOEXa6mJpqXn9WFRh87Q87vS5F5iiacJ8PTL99rT3y3Ml8AHh4pm+Y2a1mtsnMNnV3ayW6dPbigT5A4Z5qSoLZvOeyRn75XDs7O/u9LkfmIJpwn2nCdMZmWDN7J7AS+MpM33fO3e6cW+mcW1lZWRl9lZJyth3oI5jtpzhPa7inmurCXLKzfHzk7i1elyJzEE24twP10x7XAQdPPMjMrgH+D3CDc06X2zPctgPHmF+ii6mpKJiTxZollWzv6Ofp3T1elyNnKJpw3wgsMbMFZpYN3AQ8OP0AM7sQ+C6Twa7FoTPcyHiIXZ39upiawi5fVEFFQQ4fved5uvs1VktFpw1359wEcBvwCLAduNc595KZfd7Mbogc9hWgAPiFmT1vZg+e5HSSAXZ09DMRdppvT2HZWT5uXlXPseFxPnbv84S13nvKieruEufcQ8BDJzz3mWlfXxPjuiSFbYtcTNXIPbXVFufxrzeczad+uY1vP7Gb265e4nVJMgu6Q1Vi7sUDfRTnBSjVhtgpLxR2nF9XzL/9fidrtddqSlG4S8xtO9DHufOLdTE1DZgZb71gPtVFuXzwp5t5ru2I1yVJlBTuElODoxPs6OjnvLpir0uRGMkJ+Hnf6ibyAn5u+d46th865nVJEgWFu8TUpn1HmAg7hsdCXpciMVSYG+D9VywgJ8vPu76/gZbuAa9LktNQuEtMrWvpxWfQUB70uhSJsdJgNu9b3YRzjnfcvo4dHbqDNZkp3CWm1rX0UlcaJCdLe6amo6rCXO754KX4DN5x+7Nsa+/zuiQ5CYW7xMzA6ARb2/tYWJHvdSkSRxtaj/CuS5sw4JbvrdNdrElK4S4xs2nvYUJhx4JKhXu6K8vP5tY1i6guzuWd31/PFx/ezthE2OuyZBqFu8TMupbDBPxGY5nCPRMU5wV45yWNrGws47t/auGt33panTRJROEuMbOupZfz60rIztKPVabIzvLxtgvn865LG+k4NsL1/7mW//WLF7TZRxLQv0KJiYHRCbYd6OPSheVelyIeOKu2iH+8ahFXLKrggS0HeN1Xn+CLD2+nf2Tc69IylsJdYmJjZL5d4Z65gtlZvOncWj527VLOmVfEd//UwmVffIx7N+3XwmMeULhLTKxr6SXgNy5uLPW6FPFYSTCbt19czz9etYjSYIBP3reVq//tiaTdePuu9W3H/6QThbvExLo9vVxQX0JetvrbZVJdaZB/uHIRf3VRHQeODvPNx3ezeZ/WpkkUhbvMWc/AKNsO9HGZpmTkBGaTn+b+4cpFBPw+brr9We5ctw/nNE0Tbwp3mbPfbj1E2AFmaffRVmKjtjiPD1+1mNWLK/iXB17kE/dtZWRc6w/FU1SbdYicygPPH6CmKJeaolyvS5Eklpft55qzqvGZcd/mdpo7jvGdv72Y+jKtQxQPGrnLnOzrHWRL21EuqC/xuhRJAT4zrjmrmndf2si+3iGu/4+1/PDpViZCurs11hTuMie/fv4gZmj9dpmV5bVFfHDNImqKc/nX/36ZN3/jKTa0HvakluaOYzzXdoRwml0H0LSMnDHnHA88f4BVTWWUBLO9LkdSTFl+Nu+9vImXDh7j8R1d/M13n+V1yyr5n9cu5by6xHwSPHh0mLs3tDEecqzd1U19WZArl1Ym5LXjTSN3OWMvHjhGS/cgb71wvtelSIoyM86ZX8wH1yzijSuq2bL/KDd882k+8KONPLmzO+43P33x4WacgxsvmMd4yPGeH2zgb/7r2bRoDNDIXc7Yr7YcINvv4/pzavnttkNelyMpLDvLx5XLqrhkYTnP7OnlmT09PNrcRUNZkJtW1XPD+fOoK43thdcNrYf57xcOcvXyKi5ZUM7FDaX899aDbNx7mDVpMHpXuMsZmQiF+cWm/SyuKlCwS8zkBvxcvbyKNUsqKA4GuGt9G1/+3Q6+/Lsd1Jfmce78Yr7w1nMoL8iZ0+uEwo7PPfgS84pzWbNkMsiz/D6uXl7N5n1HWN/ay20sjsVb8ozCXc7I/c+10z86wUotNyBxkOX3MTga4sYL5vPaJZVsO9DHtgNHeejFDv64vYs3n1fLOy9r5KKGM/v5+/nGNl4+dIxv3Hwh/SMTx58vzgtwVm0Rm/YeYWQ8RG4gde+41py7zNrIeIiv/3EX9aV5LKsp9LocSXNl+dlcubSS2163hI+8fgkXNZbw222H+MtvP8Pbv/MMj27vnNXcfEv3AF96uJlVC8p483m1r/r+pQvLGR4P8Zutqf2JVOEus3bnun0c6hvhDWfXYGZelyMZpLoolxvOn8+nrlvOm8+rZXf3AB/48SYu+9Kj3L+5/bS7QfWPjHPrTzcTCjuuWlrJ3Rv2v+qYhRX5VBbk8NN1++L1NhJC4S6zMjA6wbef2MMViytYVFngdTmSoXICfi5fVMHHr13GX19cB8DHf/ECV37lce5Y28KRwbFX/Z1w2PGxe1+gtWeQW1Y1nLR918y4ZGEZL+w/ytb2o3F9H/GkOXeZlTvWtnB4cIxPvHEZLx3UlmriLb/PuLChlAvqS9jZOcCTu7r5v7/dzhcfbubyReVcf24tZfnZ9A2Ps3nvEf7wcieffcsKcrJOPZd+UUMpjzV38ZNn9/HVv04R7ij6AAAHn0lEQVTNu68V7hK1jr4R7ljbynVn13B+fYnCXZKGmbGsppBlNYUcPDocuQDbx9pdPa847uZVDbz38qYZp2Omyw34eduF8/nF5nY+/ablc+7O8YLCXaIyNhHmQz/bzNhEmBW1RWlxk4ekp3klecwryeMNK6rp7B8lFHbkBfwEs/3kBvynDfYp71u9gJ+tb+POdW185Jolca469jTnLlH5wm9eZkvbUf7q4joqClNvFCOZx8yoKcplfkkeZfnZs25r3NB6mGXVhdy+tiUllydWuMtp3b+5nZ+u28etaxZy7nwtECaZY/XiCgZHJ3jw+YNelzJrCnc5pWf39PLP929lQUU+9TG+/Vsk2S2qzKemKJc7nmpJud2jFO5yUo+81MF7friBsvxsbl7VgN+nnnbJLGbGFYsr2Nk58KqLs8lO4S4zumdjGx+6czNnzyvi1jULKcjRtXfJTOfVFVNZmMP31qbW6F3hLq8wMDrBvzywjX++fxuLKgu48fz5BLMV7JK5svw+/v61C1i7q4efb4yu0yYZ6F+tHPf07h4+ed9WDvYNc8XiCt5wdjVZPv3+F/nAFQt5ancvn/31S6yoLeL8FNhWUv9yM1wo7Hh0eyfv+cEG/vaO9YxOhLj1tQu5/txaBbtIxD0b97NmcQXBHD8funMzh2dY3iDZRDVyN7PrgP8A/MAdzrkvnfD9HOAnwMVAL/AO59ze2JYqsTI8FmLD3sM8s7uHh1/soO3wENVFOVy7oporFlcQ8CvURU4UzMniby9p5HtrW/jQnZv5z5svpLoo1+uyTspOd4HAzPzATuBaoB3YCNzsnHt52jH/CJznnPsHM7sJeJtz7h2nOu/KlSvdpk2b5lq/nMLIeIgDR4fZf3iIlu5BmjuO0dzRT/OhfsZCYfxmNFUEWbWgnBW1ReqGEYlCXraPT963lSyfj1vXLOTWNQvJT2DDgZltds6tPO1xUYT7ZcDnnHNvjDz+NIBz7ovTjnkkcsyzZpYFdACV7hQnz+Rwd84RdpP/dYBzEHaOsHOEwo7RiTAj4yFGJ8KMjocZnYh8PRFmNPJ8/8gEfcPj9A2Pc2wk8t/In6nnjw6PM/3/QDDbT21xLvNK8lhUWUBTeT7ZWRqli8xW78Aoj7zcyYsH+igJBljVVMaFDaWcX1dMVVEOJcFsSvICZMXhU3C04R7Nr5v5wPRLxO3AJSc7xjk3YWZ9QDkQ88bQ7z/Vyr//fsfxx1PhOPn1yX9RGacflU4tTR5tt5NjMqRxf/56+u+zE9c6nwr1WPKbkZvtJy/gIy/gJy/bT1FegOqiXApysygLZlMazKa8IJuCnCytvy4SA+UFOdyyqoG23kHWtx5mV9cAv3+581XHmUHA58PvM6b/0/vsW1bwjtc0xLXGaMJ9pjQ4MaKiOQYzuxW4NfJwwMx2nHjMaVQQh18YSU7vOXNk4vvOyPd80xfouenM/35jNAdFE+7tQP20x3XAiQstTB3THpmWKQYOn3gi59ztwO3RFDYTM9sUzceRdKL3nDky8X3rPcdPNBNCG4ElZrbAzLKBm4AHTzjmQeA9ka/fDjx2qvl2ERGJr9OO3CNz6LcBjzDZCvkD59xLZvZ5YJNz7kHg+8BPzWw3kyP2OXziEBGRuYqqf8c59xDw0AnPfWba1yPAX8e2tBmd8ZROCtN7zhyZ+L71nuPktK2QIiKSetTkLCKShlIi3M3sB2bWZWYvel1LophZvZk9bmbbzewlM/uI1zXFm5nlmtkGM3sh8p7/1euaEsXM/Ga2xcx+43UtiWBme81sm5k9b2YZcTejmZWY2X1m1hz5d31ZXF8vFaZlzGwNMAD8xDl3jtf1JIKZ1QK1zrnnzKwQ2Ay8dfqyD+nGJu+wynfODZhZAHgK+Ihzbp3HpcWdmX0MWAkUOefe7HU98WZme4GVzrmM6XE3sx8Da51zd0Q6D4POuaPxer2UGLk7555khr75dOacO+Scey7ydT+wnck7gdOWmzQQeRiI/En+0cccmVkd8BfAHV7XIvFhZkXAGiY7C3HOjcUz2CFFwj3TmVkTcCGw3ttK4i8yPfE80AX8wTmX9u8Z+DrwSSDsdSEJ5IDfm9nmyJ3r6W4h0A38MDL9doeZ5cfzBRXuSc7MCoD7gY865455XU+8OedCzrkLmLwTepWZpfU0nJm9Gehyzm32upYEW+2cuwh4E/DhyNRrOssCLgK+45y7EBgEPhXPF1S4J7HIvPP9wM+cc7/0up5EinxkfQK4zuNS4m01cENkDvrnwNVmdqe3JcWfc+5g5L9dwK+AVd5WFHftQPu0T6L3MRn2caNwT1KRi4vfB7Y75/7d63oSwcwqzawk8nUecA3Q7G1V8eWc+7Rzrs4518Tknd2POefe6XFZcWVm+ZEmASJTE28A0roTzjnXAew3s2WRp14PxLU5IiX2UDWzu4GrgAozawc+65z7vrdVxd1q4F3AtsgcNMD/jtwtnK5qgR9HNojxAfc65zKiNTDDVAO/iiw/nQXc5Zz7nbclJcQ/AT+LdMq0AO+L54ulRCukiIjMjqZlRETSkMJdRCQNKdxFRNKQwl1EJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQN/X/gniwh2Tt7rQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ill_mod = np.clip(ill, a_min = 0, a_max = 100)\n",
    "sns.distplot(ill_mod.flatten(),  bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 553,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.distplot(ill_mod.flatten(),  bins=100, hist_kws=dict(cumulative=True), kde_kws=dict(cumulative=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 554,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2a024a01978>"
      ]
     },
     "execution_count": 554,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8nGd16PHfmdG+74u1ed/jxImzERoSSKgDIWG7kKSF0g+Qy+eSsnaB2zaFtKWUlnJbGii+wO0CSYBQqKEOZkmchSSO5diO413etO/baBtpZs79QyMzkWVrLM/MO/PqfD8fEc3o1TtHRjp6dJ7zPI+oKsYYY9zF43QAxhhjYs+SuzHGuJAld2OMcSFL7sYY40KW3I0xxoUsuRtjjAtFldxFZKuIHBORJhH5zBwfrxeRp0Rkn4i8IiJviX2oxhhjoiXz9bmLiBc4DtwOtAJ7gHtV9XDENduAfar6dRFZD+xQ1aVxi9oYY8xFRTNyvw5oUtVTqjoJPAbcPesaBQrC7xcC7bEL0RhjzKVKi+KaGqAl4nErcP2saz4H/FxE/gDIBW6b76ZlZWW6dOnS6KI0xhgDwN69e3tVtXy+66JJ7jLHc7NrOfcC/6qqXxaRG4H/EJGNqhp6zY1E7gfuB6ivr6exsTGKlzfGGDNDRM5Gc100ZZlWoC7icS3nl10+CHwfQFVfALKAstk3UtVtqrpFVbeUl8/7i8cYY8wCRZPc9wCrRGSZiGQA9wDbZ13TDLwJQETWMZ3ce2IZqDHGmOjNm9xVNQA8AOwEjgDfV9VDIvKQiNwVvuzTwIdF5ADwKPABte0mjTHGMdHU3FHVHcCOWc89GPH+YeCm2IZmjDFmoWyFqjHGuJAld2OMcSFL7sYY40KW3I0xxoUsuRvX801MccMXfsV3Xoxq7YcxrhBVt4wxqeiR3c0A7DnTT+fwBF/5xXHeeXUNORn2bW/cz0buxvX2nh2gJDeDvtFJ/v0FG72bxcGSu3G1bt8Ezf1jfOQNy3nD6nK+8fRJRvwBp8MyJu4suRtXe/nsIB6Bt2+u4ZO3r2ZgbIp/e/6M02EZE3eW3I1rBUPKvuYB1lTmU5GfxVV1RbxpbQXbnjnF8MSU0+EZE1eW3I1rnejy4fMHuKah5Nxzn7x9NUPjUzwanmw1xq2sbcC4yiMRSXtv8wC5mWmsqco/99zGmkKWluZwoHXQifCMSRhL7saVxieDHOkY5nUryvB65DVJf2VFPie6RhyMzpj4s7KMcaX2oXFCCqsq8s772MqKPM70jTIVDM3xmca4gyV340pdwxMAVBZmnfexVRV5TAWVs31jiQ7LmISx5G5cqXNogpwML/mZ51ceV4ZH803dVpox7mXJ3bhS1/AElQVZiJx/vvuKc8ndl+iwjEkYS+7GdUKqdA37qZqjJAOQl5nGksIsG7kbV4squYvIVhE5JiJNIvKZOT7+FRHZH347LiLWZ2YcMzA6yWQwRFXB3MkdYGVlPicsuRsXm7cVUkS8wMPA7UArsEdEtofPTQVAVT8Zcf0fAJvjEKsxUZmZTL1oci/P46XTfYRCisdzfunGmFQXzcj9OqBJVU+p6iTwGHD3Ra6/F3g0FsEZsxCd4eReUZB5wWtWVeYxMRWibXA8UWEZk1DRJPcaoCXicWv4ufOISAOwDHjy8kMzZmE6h/2U5GaQmea94DWrrGPGuFw0yX2uv1n1AtfeAzyuqsE5byRyv4g0ikhjT09PtDEac0m6hqY7ZS5mph3yhHXMGJeKJrm3AnURj2uB9gtcew8XKcmo6jZV3aKqW8rLy6OP0pgoTQVD9I74L1pvByjKyaAsL9NG7sa1otlbZg+wSkSWAW1MJ/D7Zl8kImuAYuCFmEZozDwi943p9vlRuGAbZKSVFbnWMWNca96Ru6oGgAeAncAR4PuqekhEHhKRuyIuvRd4TFUvVLIxJu66hsLbDlxkMnXGqop8mrpHsG9Z40ZR7QqpqjuAHbOee3DW48/FLixjFqZzeII0j1CaO39yX1mRh28iQLfPP2+N3phUYytUjat0DU9QkZ+JN4redbd0zDyyu/k1pSljwJK7cZnO4Ymo6u0Q0THTZR0zxn0suRvXGPUH8E0Eoi6xlOdnUpCVRlNPao/cjZmLncRkXKPb5weYN7lHljAaSnNp6bdVqsZ9LLkb1+gdmU7u5XnzT6bOqC3O5vhFyjKRvwjuu75+4cEZk2BWljGu0evzk+YRCnPSo/6c2uJs2gbHrR3SuI4ld+MaPSN+SvMy8MxxQMeF1BRlMzEVom90Mo6RGZN4ltyNa/SOTFJ2CSUZgNriHABaB6zubtzFkrtxhWBI6R/1X3pyL8kGoHXADss27mLJ3bjCwNgkIb20yVSYLssAtNnI3biMJXfjCr3hNsiyvIxL+rz8rHQKs9OtLGNcx5K7cYWecBtkWf6ljdxhumPGyjLGbSy5G1foHZkkJ8NLTsalL92YTu42cjfuYouYTMqKXGDUO3Lpk6kzaotzePZEL6qKXEIbpTHJzJK7cYXeET+rKvIX9Lk1RdmMTQYZGJuiJPfSavbJxFbTmkhWljEpb2IqiG8iQPklTqbOqC22dkjjPpbcTcrrG5leXbqQyVSwhUzGnSy5m5R3rlNmgTX3mmLrdTfuY8ndpLzeET8ClC6wXl6YnU5+VpqVZYyrRJXcRWSriBwTkSYR+cwFrnmPiBwWkUMi8khswzTmwnpH/BTlpJPmXfhYpbY4x8oyxlXm7ZYRES/wMHA70ArsEZHtqno44ppVwGeBm1R1QEQq4hWwMbP1+vyUL7DePqO2OJvmPhu5G/eIZqhzHdCkqqdUdRJ4DLh71jUfBh5W1QEAVe2ObZjGzE1VF7Qb5Gy2r7txm2iSew3QEvG4NfxcpNXAahH5tYi8KCJb57qRiNwvIo0i0tjT07OwiI2JMDwRYDIYikFyz2HEH2BofCpGkRnjrGiS+1xL9mYPb9KAVcAtwL3AN0Wk6LxPUt2mqltUdUt5efmlxmrMeXovs1NmxszukFZ3N24RTXJvBeoiHtcC7XNc81+qOqWqp4FjTCd7Y+Kqe3gCgIoY1NzBFjIZ94gmue8BVonIMhHJAO4Bts+65sfArQAiUsZ0meZULAM1Zi7dPj9Z6R7ysy5vJ406W8hkXGbenwhVDYjIA8BOwAt8W1UPichDQKOqbg9/7M0ichgIAn+kqn3xDNwYgK5hPxX5WZe94VdBdhr5mWmuTO6258ziFNVwR1V3ADtmPfdgxPsKfCr8ZkzcRCYqgG7fBBuWFFz2fUWEGtvX3biIrVA1KWvEH2BsMkhFflZM7ldfkkNzvyV34w6W3E3KOjeZWnB5k6kzZpK79bobN7D93E3K6g6fm3o5I/fIMk9DaQ4TUyG6fX4qC2Lz14AxTrGRu0lZXcMTZKZ5KLjMTpkZ9aW5AFaaMa5gI3eT9GZPos6YGWHH6mi8+pLpdsizfWNcu7QkJvc0xik2cjcpq3t44rIXL0WqKcrGI9DcNxqzexrjFEvuJiWN+AOMTgapiGFtPCPNw5KibCvLGFew5G5SUrcvNtsOzFZfksPZJE3uU8EQQ2NThEKv7eZRVevwMeexmrtJSheqs8/oHp7ulIl1V0tDaQ4/P9QV03vGwsDoJO/6l+c51TOKR6AoJwOvRxifDDLqD5CbmcbqynzWVOWzqiLP6XBNErDkblJSt88f006ZGfUlufSNTjLiD5CXmRw/HoFgiI89to/W/nE+fftq/IEQA2OThFTJTk/jVO8I/aOTHO4Y4uXmAbLSPXg9wu9cX39Zp1OZ1JYc373GXKKZydRYdcrMaCid7php7htjfQy2NYiFv/v5MZ490cuX3rWJ91xbd97HZ/7KCYaU5v4xnjrazV9sP8SjLzXzl2/fmOhwTZKwX+smJXX7/DGdTJ0x0w7Z3J8cHTM/OdDON54+xftuaJgzsUfyeoRlZbn8/k1L+ZffvRrfRID7/u+Ltl/OImXJ3aScMX+AEX8g5pOpAPWlv+l1v5hHdjefe4uXqWCIP/vxq1zTUMyf37k+6s8TEbZurOa/P/Z6KvKzePSlZsYng3GL0yQnS+4m5XT54jOZClCQlU5xTnpSdMw0nhlgaHyK+29eTkbapf+oFuVk8M/3bWZ4PMDje1uso2aRseRuUk7H0PSe6/Ha/6W+JIeWJEjuTx7tIsPr4fUryxZ8j831xWzdWMWRTh/PNfXGMDqT7Cy5m5TT3D9GYXY6hdnpcbl/fWnuvGWZRPjV0W5uWFFK7mV27bxuRSnrqwvYeaiTzqGJGEVnkp0ld5NyWvrHqAtPfMZDQ0kObYPjTAVDcXuN+ZzuHeVUzyhvWltx2fcSEe7YWIUq/O8fHYz7XIFJDtYKaVKKb2KKgbEpblyeHbfXqC/JIRhSOgYvfZQbqyPtnjzaDcAbY5DcAUrzMlldmc9Lp/u5ZU05aR4b17ldVP8Pi8hWETkmIk0i8pk5Pv4BEekRkf3htw/FPlRjoKV/ut4ez5H7uY4ZB9shnzzaxerKPJ490RuzkfaNK0oZ8Qc41DYcgwhNsps3uYuIF3gYuANYD9wrInP1ZX1PVa8Kv30zxnEaA0DLwBheEZYUxW/k3hBlO2S8+Cam2H2qnzeurYzpfVdW5FGam8ELp+zs+sUgmpH7dUCTqp5S1UngMeDu+IZlzNya+8eoLsoiPY7L6ivzs8hI8zi2O+SzJ3oJhJTb1sWmJDPDI8INy0tp7h+jbXA8pvc2ySean5AaoCXicWv4udneJSKviMjjIjLnUjoRuV9EGkWksaenZwHhmsUsGFJaB+I7mQrg8Qh1xdk0OzRy/9WRbopy0tlcXxzze19dX0y6V3jRRu+uF01yn2vzjtmrIX4CLFXVTcAvgX+b60aquk1Vt6jqlvLy8kuL1Cx6XcMTTAWV+uL4JneA5eV5nOj2xf11ZguGlJ+92sHS0ly+t6dl/k+4RNkZXjbXFXOgZZDhiamY398kj2iSeysQORKvBdojL1DVPlX1hx/+X+Ca2IRnzG+0hPdIiffIHWDDkgJO9Y4yGUhsO+T+lgFGJ4OsrcqP22tcWVdEIKQ832SjdzeLJrnvAVaJyDIRyQDuAbZHXiAi1REP7wKOxC5EY6Y1942Rm5lGcU58Fi9F2rikENXfrIZNlF8c7sYrwurK+CX3+pIcMtI8PHvCSqNuNm+fu6oGROQBYCfgBb6tqodE5CGgUVW3Ax8TkbuAANAPfCCOMZtFqmVgjPqSnJhv8ztjpt3wvuvr2VhTCED74DgNpblxeb25/PJIF8vKc8lK98btNbweYUVZLs+esO0I3CyqRUyqugPYMeu5ByPe/yzw2diGZsxvjPkD9I5Mck0cJhnnUlmQSWluBu0LWMi0UKd7R2nqHuFtm6rnv/gyrazM5ycH2jnbN5rQX14mcWyFqkkJLQPxX7wUSUTYUFPIia7ETKo+srv5XJlkbfXch4TEavUrcO4ovmdO9PI+S+6uZGuQTUo42TOC1yPUJqBTZsbGJQV0DU8QSNAeM0c6fFQXZlGckxH31yrNzaA4J51HbZ8Z17LkblLCsS4fy8pyF7Sv+UJtrCkkpNA17J//4nnMt4XAmD/A2b5R1lYl5mg/EWFlRT4ne0YIhmyfdzeysoxxVDSlhv7RSXp8fq5bWpKosIDpjhmYnlStKY7NdgcXSvDHunwosK46fl0ys62qyGPPmf6k2LvexJ6N3E3SOxaue6+JY+/3XOpKsslK99CWgHbIIx3DFGSlxXXPnNlWlOchwInukYS9pkkcG7mbpHe800dpbgZlebE/M/ViRIQlhdm0x3kflslAiOPdI1xVW4QnyjbPWNTJszO81BZn0+TASlwTfzZyN0ltMhDiZM9IwkftM5YUZdM5NBHXuvRLp/uYDITYXF8Ut9e4kFWV+bQOjDM4Npnw1zbxZcndJLXTvSMEQsqaOK7YvJglRVkEQkrPyOVPqs5lMhDi6RO9rCzPc6TffGlpLgocbBtK+Gub+LLkbpLa0U4f6V5hWZkzvdhLCqdr4PEqzbx0uo9RfyBmJy5dqurC6UPGj3ZYacZtLLmbpKWqHO/ysbI8j7Q47t9+MWX5maR7JS7JPXLUvtShX165mWkUZKVxpMNOZ3IbS+4maXX7/AyMTbHaoXo7TB9wUV2YHZfDLXaHR+1vivGhHJeqqjCLI502cncbS+4maR2faYF0qN4+Y1lZLi39Y/SPxm7ScXwyyDMnellZ4UytPVJVwXTHTKK3NzbxZcndJK2TPSOU52VSlIDl+BdzRXil6s9e7YzJ/QKhEN/dfZaJySBvXh/bc1IXoqowi6mgcqrX+t3dxJK7SUqBUIjTvaOsqHB+U6vqwizK8jL4yYH2+S+eh6ry431tnOod5Z1X1yR0r5wLsUlVd7LkbpJSa/84U0FlRXme06EgIlxRU8Tu0310+y5vC+CnjnXzcvMgb1pbEZczUheiLC+TDK/HJlVdxpK7SUone0YQYHmZ88kdYFPtdGnmiYMLK80EQ8rOQ5388kg3m+uKHGt9nIvXI6ysyDs3qTrfJmcmNVhyN0npZM8IS4qyyc6I34lEl6KyIIs14QMuLlXviJ//9/xpnj7ew7VLS3jH1TVxO01qodZVF3DURu6uYsndJJ2xyQAt/eOsKHe+3h7pzk3VNJ4duKSe92eO9/C2rz5Hc98Y77q6lndsriHNk3w/duuq8+n2+emL00pck3hRfZeJyFYROSYiTSLymYtc924RURHZErsQzWKz58wAQU2OenukO69cAsCOgx3zXuubmOLjj+3j/d9+iex0Lx95wwquaUiOGvtc1oVPfzpq/e6uMW9yFxEv8DBwB7AeuFdE1s9xXT7wMWB3rIM0i8vzTb14RRzv/55tWVkuG5YU8MOX2y66kdjRjmG+8svjPHGwk0/ctoonPvFbCd3KdyHWhheK2aSqe0Qzcr8OaFLVU6o6CTwG3D3HdX8JfAlI3InCxpWeP9lHXUlOQk9dijQzoTjXpOKHfmsZRzqG2Xlo7onVvWf7+c7us5TkZvDEJ36LT9y2msy05Jg3uJjSvEwq8jM5Yu2QrhHNT08N0BLxuDX83DkishmoU9WfXuxGInK/iDSKSGNPT88lB2vcb3Bsklfbh5Kiv30u79hcy/tvbOC5pl72NQ+ce15V+dquJn74chvLy/P48OuXJ11ZaT5rqws42mkjd7eI5rCOuab1z/1NKiIe4CvAB+a7kapuA7YBbNmyxQ5uNOd58VQfqrAiSVog5/Lnd67n2RO9/GhfG5lpXv7l6ZP89JV2Xm0bZlNtIe++pjYpJ00v5pHdzagqRzt9BEOK15Nc3Tzm0kWT3FuBuojHtUBkP1g+sBHYFW7vqgK2i8hdqtoYq0DN4vDiqX6y073UliRvjTrd6+G+6+r52q4mvrP7LABX1hXx+bs24PVI1KcpJZvqwiyC4b3rqwqynA7HXKZokvseYJWILAPagHuA+2Y+qKpDQNnMYxHZBfyhJXazEAdaB7mipjDpR765mWn8/k3LON7l40+2rqWuZHobgWRZ/LOQOCrDCb3HZ8ndDeZN7qoaEJEHgJ2AF/i2qh4SkYeARlXdHu8gjbtcKPFMBUMcbh/mfTc0JDiihSnLy6QsL/NcYk91pbmZCNDjmwAKnQ7HXKaoDshW1R3AjlnPPXiBa2+5/LDMYnS8y4c/EOKK2kJG/UGnw1l0MtI8FOak0zti56m6QXL/7WsWlYOt0+d4Xlmb+IOizbTyvEx6fLZK1Q0suZukcaB1iIKsNBpK3VHmSEXl+Zn0jPhRtWa2VGfJ3SSNg22DbKotSrpNtRaTsrxMJgMhhicCTodiLpMld5MUpoIhjnb4uKLWJvKcVJ6fCWClGRew5G6SQufQBIGQsqnGkruTziV32x0y5VlyN0mhNbyN7qY6m0x1Un5mGplpHnpt5J7yLLmbpNA2ME5pbgZLCm3xjJNE5NykqkltltxNUmgdGGNTbaFNpiYBa4d0B0vuxnH+QJAen58rrL89KZTlZzI0PsXYpHXMpDJL7sZx7YMTKHCldcokhfK86UnVUz2jDkdiLocld+O4tvBkqrVBJoeZjpmTPSMOR2IuhyV347i2gTEKstKoyLfJ1GRQmpuBYCP3VGfJ3TiufWgi6c8YXUzSvB6KczNs5J7iLLkbR00GQvT6/Jbck0x5XiYnbeSe0iy5G0d1Dk9Pplp/e3Ipz8/kdO8IoZBtIJaqLLkbR7WHJ1OrbeSeVMrzMpmYCtE+NO50KGaBLLkbR7UPjpOd7qUoO93pUEyEsnDHTFO31d1TlSV346iOoQmWFGXZytQkU1kwndyPdfocjsQsVFTJXUS2isgxEWkSkc/M8fGPiMhBEdkvIs+JyPrYh2rcJhhSOocnWFJoJZlkk5ORRnVhFkc6hp0OxSzQvMldRLzAw8AdwHrg3jmS9yOqeoWqXgV8CfiHmEdqXKfbN0EwpFZvT1Lrqws40mEj91QVzcj9OqBJVU+p6iTwGHB35AWqGvnrPRewKXYzr47BCcA6ZZLVuuoCmnpGmJiyw8pTUTTJvQZoiXjcGn7uNUTkoyJykumR+8fmupGI3C8ijSLS2NPTs5B4jYu0D42T7pVzk3cmuayrLiAYUptUTVHRJPe5ZrrOG5mr6sOqugL4E+DP5rqRqm5T1S2quqW8vPzSIjWu0z44QVVBFh6bTE1K65cUAHC43eruqSia5N4K1EU8rgXaL3L9Y8DbLyco434hVTqGxm1lahJrKMkhJ8PLYZtUTUnRJPc9wCoRWSYiGcA9wPbIC0RkVcTDtwInYheicaOB0Un8gZB1yiQxj0dYU5VvHTMpKm2+C1Q1ICIPADsBL/BtVT0kIg8Bjaq6HXhARG4DpoAB4PfiGbRJfe1D05Op1UWpM5n6yO5mp0NIuPXVBWw/0I6q2lqEFDNvcgdQ1R3AjlnPPRjx/sdjHJdxuY7BcTwClQWpk9wXo3XVBXx3dzNtg+PUFuc4HY65BLZC1TiiY2iCivws0r32LZjM1lVPT6pav3vqiWrkbkysdQyNs6I8z+kwzDzWVuUjMt0xc/v6SqfDSVqRJbv7rq93MJLfsGGTSbgxf4DhiQBVtngp6eVmprG0NNcmVVOQJXeTcB3D05OpltxTw7rqfI50WnJPNZbcTcJ1zHTKWBtkSlhXVcDZvjF8E1NOh2IugSV3k3CdQxPkZ6aRl2lTPqlgZqWqbf+bWuyny8TVzERT5CRT59D4oijJuKUvfsOSQgD2NQ+yZWmJw9GYaFlyNwk1FQzR5fPzugr3dMq4JYlfSFVhFmsq8/nV0S4+fPNyp8MxUbKyjEmo072j03u4L4KRu5vctr6CPWcGGBybdDoUEyVL7iahZlrqqmwyNaWEQtMnZ31hxxGnQzFRsuRuEupwxzBej1CeZ3u4p5Ka4mzyMtNspWoKseRuEupIh4+K/Ey8HtuEKpV4RFhblc/xLh+TgZDT4ZgoWHI3CXW0Y9jq7SlqXXUB/kCI3af7nA7FRMG6ZUzC9I346fb5rZ0uRczuAlpRnke6V/jl4S5+a5WdpJbsbORuEmamXltl2/ympIw0DyvL8/jlkW5Uzztp0yQZS+4mYWY6Zawsk7rWVRfQNjjOUVutmvQsuZuEOdI5TGVBJrm27UDKWlOVj0fgP19udToUMw9L7iZhDrcPs7aqwOkwzGXIz0rnytoi/t+vz/D1XSedDsdcRFTJXUS2isgxEWkSkc/M8fFPichhEXlFRH4lIg2xD9WkMn8gSFP3CBtrLLmnutvWVaIKTx7tdjoUcxHzJncR8QIPA3cA64F7RWT9rMv2AVtUdRPwOPClWAdqUtvxzhECIT23CZVJXcW5GVy7rIS9Z/s53TvqdDjmAqIZuV8HNKnqKVWdBB4D7o68QFWfUtWx8MMXgdrYhmlS3aH2IQA2LLGRuxvcuqacNI+HL//8mNOhmAuIJrnXAC0Rj1vDz13IB4En5vqAiNwvIo0i0tjT0xN9lCblvdo+RH5mGnXFOU6HYmIgPyud160s5aevdPBq25DT4Zg5RJPc51onPmeTq4j8LrAF+Lu5Pq6q21R1i6puKS+3RRCLyaH2YdYvKcBj2w64xs2ryinOSeez/3nQtiRIQtEk91agLuJxLdA++yIRuQ34U+AuVfXHJjzjBiFVjnb4rN7uMlnpXv72XZs42DbEl3521OlwzCzRJPc9wCoRWSYiGcA9wPbIC0RkM/ANphO7TaGb1+j1+RmfClq93YXevKGK99/YwDefO81T1j2TVOZN7qoaAB4AdgJHgO+r6iEReUhE7gpf9ndAHvADEdkvItsvcDuzCLUPjQOwwdogXeeR3c2sKM+jqiCLT//gAF3DE06HZMKi6nNX1R2qulpVV6jqX4efe1BVt4ffv01VK1X1qvDbXRe/o1lM2gcnyAzvS2LcJ93r4Z5r6xifDPKxR/cRCFr9PRnYClUTd+1D46ytyifNa99ublVRkMUX3rmR3af7+dJOa49MBvbTZuJKVekYnGC9Taa63js21/K+GxrY9swpdhzscDqcRc+Su4mrwbEpxqeCtu3AIvHnd65nc30Rf/SDAzR1jzgdzqJmyd3E1bnJVBu5u94ju5t5fG8rb15fRVa6l498Zy+j/oDTYS1atveqiav2wQk8Amur8p0OxSRIYXY6/3TvZn73W7v50x8d5NqlJYhML1677/p6h6NbPGzkbuKqY2icsrxMstK9TodiEuimlWV86rbV/Hh/Oy+d6Xc6nEXJkruJG1WldWCcmqJsp0MxDvjorSu5ZU05P32lg9aBsfk/wcSUJXcTN2f7xhjxB2gozXU6FOMAj0f4ynuuIi8zje/ubsY3MeV0SIuK1dxNTDyyu/nc+zN11Zk/xxtKbSfIxSby++F9NzTwjWdO8h8vnuX9Ny4lO8NKdIlgI3cTN41n+snJ8FKRn+l0KMZBS4qyee+WetoGxvn0D/YTCs25qayJMUvuJm72nBmgoSTnXKeEWbzWLylg68Yqdhzs5G9/dhTV3yT4R3Y3n3szsWNlGRMXPT4/p3tHuWNjldOhmCTx+pVllOZl8I1nTjE2GeRzd23Aa/v7x40ldxMXe8/O1NttMtVMExH+8u6N5Gam8Y2nT9E36ucf3nOV02G5liUacgUSAAANvklEQVR3ExcvnR4gK93DkqIsp0MxSURE+Owd6yjPy+Sv/vsIXcO7uWlFGVWF9n0Sa1ZzN3HReLafq+qKSPPYt5g5X05GGu+9to7D7cN89ckT/Nf+NtuqIMZs5G5ibtQf4FD7MB+9ZYXToZgkdmVtEasq8vjVkW52n+7jQOsgQVU+8LqltqI5BmxYZWJuX/MgwZCyZWmJ06GYJJeTkcbbrlzCH7xxFQ0luXzxiaPc+ve7eHxv62s6alLBZCDEy2cHkqbV00buJuZeOtOPR+DqhmJaB8adDsckkQu1O1YWZPF7r1vKsrJcvvjEEf7wBwf4+q4m3rG5lpLcjJTYcOypY908fbyHd2+p5YblpU6HE93IXUS2isgxEWkSkc/M8fGbReRlEQmIyLtjH6ZJJY1n+lm/pIC8TBs7mEtz44pSfvzRm/jCO66gdWCcf/zVcZ4/2Zs0o+ELGfEHeP5kLwAHWgYdjmbavMldRLzAw8AdwHrgXhFZP+uyZuADwCOxDtCkFn8gyN6zA1y31PmRi0lNIsJ919fz8TetYllZLj99pYP3bnuBUz3Je/jHM8d7CASV7HQvr7QOOR0OEF1Z5jqgSVVPAYjIY8DdwOGZC1T1TPhjdjLuInes04c/EGKrLV4yCxBZtinKyeD3blzKvuZBfn64k63/+Cyfun01H3r9sqQ6j7d7eIIXT/Wxub6IyaByoDVFRu5ADdAS8bg1/Jwx5znYNkRFfiZbGoqdDsW4gIhwdUMx/+vWlawsz+OLTxzlrn/+ddKUPgC+tuskIVVuXVNBbVE2rQPj9I34nQ4rquQ+1/rgBRXAROR+EWkUkcaenp6F3MIkMf9UkGOdPt5yRTUeW1ZuYqggK53fub6e+66rp2/Uz9u/9ms+t/2Q49sIdwyN88juZq6uL6Y0L5PakumzC5KhNBNNcm8F6iIe1wLtC3kxVd2mqltUdUt5eflCbmGS2JFOH4GQcuemaqdDMS4kImysKeR/3ryC65eV8G/Pn+GmLz7Jn/7ooGMxfW9PC1OhELeuqQCgpjAbEZKiNBNNzX0PsEpElgFtwD3AfXGNyqSE2W1tB9uGKMhK4+p6K8mY+MlK93LXlTVcVVfMj/e18d3dzXT7/Dx09waqCxN76tfOQ11saSimODcDgMx0LyvL85KibDRvclfVgIg8AOwEvMC3VfWQiDwENKrqdhG5FvgRUAy8TUQ+r6ob4hq5SSoTU0GOd/m4YVmJlWRMQtSX5PDRW1fy66Zedh3v5rYvP80f/vYa0r0ePAk4kPts3yhHOob5s7eue83zV9YV8dTRblTV0e2uo5pyVtUdqrpaVVeo6l+Hn3tQVbeH39+jqrWqmquqpZbYF58jHcMEQ8oVNYVOh2IWEa9HuHl1Ob/45BvYsrSEz//kMF/fdZKW/vif2brzUCcAv73htZ1hV9YW0jc6Sdugswv4kqefyKS0g21DFGanU1tiR+qZxKsryeFff/9avnrvZobHp/j60yd5bE9zXA/m/tmrnWxYUkDdrO/5TbVFABxocXZS1ZK7uWxD41Oc6BrhiprCc38OG5NIj+xu5tGXWvBNBPjU7au5ZU05h9uHeeOXn+aPfnAg5qtcu4cneLl5kK0bzl/PsbY6nwyvh1ccnlS19eHmsu061g3AjUmwn4Yxmele3ry+iuuWltAyMMb2/e38YG8rhdnpbKotZFNNEX/426svqx6+83AXwJyL9TLTvKyrzme/w5OqNnI3l2VgbJLGMwNcs/Q3HQPGJIOinAz+5p2baPyz2/mnezdTVZDFr5t6eXhXE7f8/S7+6VcnFrzYaOernSwvy2VlRd6cH99UW8SrbUMEHdwTx0bu5rI8dbQbEc71+RqTbLIzvNx15RJGJgKMTQY43D5Ml2+Cf/jFcR5+qol3bK7h/TcuZf2SgqjuNzg2yYun+vjwzcsvOPq/sq6I/3jxLCd7RlhdmR/LLydqltzNgvWN+Hm5eYDrl5dSmJ3udDjGnGf2WoycjLRz5wxc21DC8yf7eHxvK4/taaGqIIvfv2kpb91UTW3x3I0BgWCIr+06SSCkc9bbZ1xVNz2pur950JK7ST1PHu3G6xFuWW2rjU3qqSjI4u2ba3jz+kpeaRtiX/MAf/PEUf7miaMsK8ulPD+TZWW53H/zcpaW5nKs08ef/PAVDrYNceemajbVXrjtd3lZLkU56ew9O8B7rq274HXxZMndLMi+5gH2twxy08oy8rNs1G5SV05mGjcsL+WG5aW8bkUpTx7t5tdNvTzb1MtLp/v53p4WstO9TAVDFOWk8/B9V/OWK6ouOiHr8QhX1xezt3kggV/Ja1lyN/OK/NP2vuvrGZsM8Mnv7acwO503rrVau3GP50/2kZXu5U3rKrllTQVdwxN0Dk3QMTSOiPDVezefaxy40KlSM65pKObJo90Mjk1SlJP4ZgNL7uaSfWHHEc72j/HBm5bZQcbGtbweYUlRNkuKspneWQWeeLUz6s+f2WNpX/MgtzowCLJWSHNJdh3r5jsvNvOh1y9jefncbWDGGLiyrhCvR9h71pnSjCV3E7WBsUn++PFXWFOZz6ffvMbpcIxJajkZaayvLrDkbpJb1/AE33j6JBNTQf7PPVdZOcaYKFzTUMz+lkECwcSfQGrJ3cyruW+Ubc+cQoHvf+RG1lVHt9jDmMVuc30R41NBjnb6Ev7altzNRW0/0M63fn2anAwvH7l5BWurLLEbE61rwmcJv+xAS6R1y5hzIlu73n1NLX/934f5txfO0lCSw33X11s/uzGXqKYom8qCTPaeHeD9Ny5N6Gtbcjevoaqc7RvjvdteYF/zIB96/TIaSnPx2ulKxlwyEeGahmJHJlUtuRtgek/2F0718dLpPrqG/RRkpfG137mat1xRPe9iDWPMhV1dX8yOg510DU9QWZCVsNeNKrmLyFbgH5k+Q/WbqvrFWR/PBP4duAboA96rqmdiG6q5XLNXmo5PBmk8288PGlvZeagTfyBETVE279xcw1+9YyM5Gfa735jLdW14o7Kv7zrJX7xtfcLOVZ33p1dEvMDDwO1AK7BHRLar6uGIyz4IDKjqShG5B/hb4L3xCNjMLzKJv+uaGgbHpmgdGOel0/10Do/TPeznq0+eoGNoAoDC7HTee20deZlp53bDs8RuTGxsqi3k/Tc28K/Pn0FV+Yu3bUjIIfLR/ARfBzSp6ikAEXkMuBuITO53A58Lv/848M8iIqrq3E71SS4UUqZCISYDIaaCr/1nUlVCOv3fofEp+kcnGRibZGwyiD8w/TkhVTwiiMDA6BQtA2O09I8xODZFz4ifQDCEPxDif//o4GvunZnmoSI/kxtXlLK0NJc1Vfm8YXU5WeleK78YEwciwufv2kCG18M3nzvNZDDEX7/9irgn+GiSew3QEvG4Fbj+QteoakBEhoBSoDcWQUb61nOn+fLPj517rAqKhv9L+H8iCAgw85eQIMz1V9HMr6GZe80l8h7nPm/mtRVCqgT1/M8XAY8IHoFQ+LpY/9qrLMikrjiHpWU5ZKV7SPN6yEjz8LrlpRTnZlBZkMXxTh9FOemICPddXx/bAIwxFyQi/Olb15GZ7uHhp07SUJrLR96wIr6vOd/gWkT+B/Dbqvqh8OP3Adep6h9EXHMofE1r+PHJ8DV9s+51P3B/+OEa4BiXr4w4/BJJIYv96wf7N1jsXz8srn+DBlWd9xCFaEburUDkbvO1QPsFrmkVkTSgEOiffSNV3QZsi+I1oyYijaq6JZb3TCWL/esH+zdY7F8/2L/BXKJZoboHWCUiy0QkA7gH2D7rmu3A74XffzfwpNXbjTHGOfOO3MM19AeAnUy3Qn5bVQ+JyENAo6puB74F/IeINDE9Yr8nnkEbY4y5uKj63VR1B7Bj1nMPRrw/AfyP2IYWtZiWeVLQYv/6wf4NFvvXD/ZvcJ55J1SNMcakHtsV0hhjXChlk7uIfFtEukXkVadjcYKI1InIUyJyREQOicjHnY4pkUQkS0ReEpED4a//807H5BQR8YrIPhH5qdOxJJqInBGRgyKyX0QanY4nmaRsWUZEbgZGgH9X1Y1Ox5NoIlINVKvqyyKSD+wF3j5rWwjXkukNOnJVdURE0oHngI+r6osOh5ZwIvIpYAtQoKp3Oh1PIonIGWCLqi6WHveopezIXVWfYY5e+sVCVTtU9eXw+z7gCNMrhRcFnTYSfpgefkvNkcplEJFa4K3AN52OxSSXlE3u5jdEZCmwGdjtbCSJFS5H7Ae6gV+o6qL6+sP+D/DHQOIP6UwOCvxcRPaGV8CbMEvuKU5E8oAfAp9Q1WGn40kkVQ2q6lVMr5q+TkQWVXlORO4EulV1r9OxOOgmVb0auAP4aLhca7DkntLCteYfAt9V1f90Oh6nqOogsAvY6nAoiXYTcFe47vwY8EYR+Y6zISWWqraH/9sN/IjpXWwNltxTVnhC8VvAEVX9B6fjSTQRKReRovD72cBtwFFno0osVf2sqtaq6lKmV4U/qaq/63BYCSMiueFmAkQkF3gzsCi75+aSssldRB4FXgDWiEiriHzQ6ZgS7CbgfUyP1vaH397idFAJVA08JSKvML3/0S9UddG1Ai5ylcBzInIAeAn4b1X9mcMxJY2UbYU0xhhzYSk7cjfGGHNhltyNMcaFLLkbY4wLWXI3xhgXsuRujDEuZMndGGNcyJK7Mca4kCV3Y4xxof8PWvnYWIfEwkIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(s_ori.flatten(),  bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 564,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 4\n",
    "\n",
    "\n",
    "L = np.exp(ill_mod)\n",
    "\n",
    "S = np.exp(s_ori)\n",
    "\n",
    "R = np.exp(ref)\n",
    "\n",
    "W = np.ones(S.shape) * 255\n",
    "\n",
    "L_ = W * np.power((L/W), 1/gamma)\n",
    "\n",
    "# S_ = (L_/L) * S\n",
    "\n",
    "S_ = L_ * R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 565,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x2a025e36cc0>"
      ]
     },
     "execution_count": 565,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0ZGd95//3t6pU2vel1VIv6t3ubu/tBWMbCJtNgIaMSQxkYCYmzmTwbybDZM7AJOEwDJzzc+aX8fwyEAiLgw1xMGYxbTBuvGHjrbvVdtu9b2p1a19a+y5VPfNHXQlZraWkKlWVpM/rnDoq3bp161ulK330PPc+zzXnHCIiIr5kFyAiIqlBgSAiIoACQUREPAoEEREBFAgiIuJRIIiICKBAEBERjwJBREQABYKIiHgCyS5gPkpKSlxVVVWyyxARWVIOHjzY7pwrnWu9JRUIVVVVVFdXJ7sMEZElxczOR7OeuoxERARQIIiIiEeBICIigAJBREQ8CgQREQEUCCIi4lEgiIgIoEAQERGPAkFERIAlNlJZVq6H912YuP+JG9clsRKR5UstBBERARQIIiLiUSCIiAigQBAREY8CQUREAAWCiIh4FAgiIgIoEGQJaeoeJBR2yS5DZNnSwDRZEroGRvjas2e4bWspfp9NLNcgNZH4UQtBloSatn4c8GrNRQZHQskuR2RZUiDIklDT3k+a3xgeC7Pv3MVklyOyLCkQZEk4197HlrJctq7K4aUz7YyMhZNdksiyo0CQlFffOUDnwCgbS7N5x9Yy+kdCVJ/vSHZZIsuOAkFS3r6ayB//DSXZbCjJZn1RFi+ebtcZRyJxprOMJGWNT3n9k4P1ZKb5WZWXAcA7t5Xy4CvnebO+C1ifxApFlhe1ECTl1bT3saEkG59FTjfduiqXjDQfdZ2DSa5MZHlRIEhK6xwYmTh+MM7MKMwK0tk/ksTKRJYfBYKktHPt/UDk+MFkRdlBOhQIInGlQJCUdq6t/y3HD8YVZQXpHBghrAPLInGjQJCUNvX4wbjC7CBjYUdb33CSKhNZfhQIkrK6B0fpHBi9pLsIIl1GABc6BhJdlsiypUCQlNXu/fdfnp9xyWNFWZFAqFMgiMSNAkFS1vhZRIXeH//JCrLSMNRCEIknBYKkrM6BEQzIz0y75LGA30deZpoCQSSOFAiSsjoHRsnPSnvL9Q8mK8wKUt+hwWki8RJVIJjZ7WZ20szOmNnnp3k83cwe8R7fZ2ZV3vL3mtlBMzvsff29Sc+5zlt+xsz+3sym/62XFauzf2Ta7qJxRdlqIYjE05yBYGZ+4OvAHcB24ONmtn3KancDnc65zcD9wH3e8nbgQ865K4BPA9+f9JxvAPcAW7zb7TG8D1mGOgdmD4TC7CAtvUMMjeqCOSLxEE0L4QbgjHOuxjk3AvwQ2D1lnd3Ag979HwPvNjNzzr3unGv0lh8FMrzWxGogzzn3inPOAQ8BH4n53ciyMTwWomdojMLsS48fjCvKCuIcNHSp20gkHqIJhEqgbtL39d6yaddxzo0B3UDxlHX+FfC6c27YW79+jm0CYGb3mFm1mVW3tbVFUa4sBw3exHWzdxlpLIJIPEUTCNP17U+dL2DWdcxsB5FupD+bxzYjC537lnNul3NuV2lpaRTlynJQF0UgFHqBUK9AEImLaAKhHlg76fs1QONM65hZAMgHOrzv1wA/Az7lnDs7af01c2xTVrD6zsgf+cKsmbuMctMDBHzG3qMtE9dOEJGFiyYQDgBbzGyDmQWBu4A9U9bZQ+SgMcCdwLPOOWdmBcAvgS84514aX9k51wT0mtlN3tlFnwJ+HuN7kWWkrmMQvxl504xBGGdmFGrWU5G4mTMQvGMC9wJ7gePAj5xzR83sy2b2YW+17wLFZnYG+BwwfmrqvcBm4G/M7JB3K/Me+3PgO8AZ4Czwq3i9KVn66jsHyM9Ku2RSu6nGZz0VkdhFdQlN59wTwBNTln1x0v0h4GPTPO8rwFdm2GY1sHM+xcrKUdc5ODFf0WyKsoPUXuwncrKaiMRCI5UlJTV0Dsx6yum4wuwgw2NhBkc0FkEkVgoESTkDI2O0980+KG3ceCuiQ91GIjFTIEjKiWYMwrjxsQg6sCwSOwWCpJy6KE45HTferdSpQBCJmQJBUk6dN4Pp+MCz2aQH/GQF/XQOji52WSLLngJBUk595wDpAR856VGdBEdeRhq9CgSRmCkQJOXUdQyypjCTaGdEz80I0Ds8tshViSx/CgRJOfVdA6wtyop6/byMNHrUQhCJmQJBUs54CyFauZkB+obHCIU1OE0kFgoESSkPvHiO7sFR2nujP2soNyONsNOppyKxUiBIShmflyiaM4zG5WVEDj639AwtSk0iK4UCQVJK10DkWEA0YxDG5WVE1m3tVSCIxEKBICllvIVQEMUo5XG5XguhtWd4UWoSWSkUCJJSugZGCfiM7KA/6ufkTHQZKRBEYqFAkJTSNTBCQVZa1GMQAAI+H9lBPy3qMhKJiQJBUkrX4Oi8uovG5WakqctIJEYKBEkpXQOjFMxy2cyZ5GUGdFBZJEYKBEkZQ6Mh+obHFtxC0GmnIrFRIEjKaOyKzHJaMI9TTsflZgRo7xvRaGWRGCgQJGU0xBAIeRlphMKOi/06jiCyUAoESRkTV0rLnH+XUZ7GIojETIEgKaOxaxAD8hZwUDlXo5VFYqZAkJRR3zVIXmYafl/0YxDG5WpwmkjMFAiSMho6Bxd0yin8roWgM41EFk6BICmjoWtwQQeUAfzedBcvnbkY56pEVg4FgqSEUNjR3D20oDEI4/Iy0+gd0pXTRBZKgSApoaVniLGwW3ALAbxrKw/p2soiC6VAkJQwMShtAaecjsvNSKNHLQSRBVMgSEqIZVDauLyMAH1DurayyEIpECQl1I8PSovhGEJuRhoOuNinU09FFkKBICmhoWuQwqw0goGF75J5GosgEhMFgqSEhs5BKgszY9qGRiuLxEaBICmhsWuQivzYAmF8ygu1EEQWRoEgSeeco6Er9hZCTnoAQ6OVRRZKgSBJ1zUwysBIiMqC2ALB7zOy0nXlNJGFiioQzOx2MztpZmfM7PPTPJ5uZo94j+8zsypvebGZPWdmfWb2tSnP+Y23zUPerSweb0iWnvFTTtfE2EIAyM8I0NytQBBZiDkDwcz8wNeBO4DtwMfNbPuU1e4GOp1zm4H7gfu85UPA3wB/OcPmP+mcu9q7tS7kDcjSV985AEBlQVbM28rPTKNJgSCyING0EG4AzjjnapxzI8APgd1T1tkNPOjd/zHwbjMz51y/c+5FIsEgMq2fvtYAwIHajpi3lZ+VNtHiEJH5CUSxTiVQN+n7euDGmdZxzo2ZWTdQDLTPse1/MrMQ8BPgK845DTFdgdr7hsnNCJCR5o95W/mZQXqHxugbHiMnPZrdOz4e3ndh4v4nblyXsNcViadoWgjTXa1k6h/uaNaZ6pPOuSuAW73bv572xc3uMbNqM6tua2ubs1hZetp6hynNSY/LtvK9U0+b1EoQmbdoAqEeWDvp+zVA40zrmFkAyAdmbf875xq8r73Aw0S6pqZb71vOuV3OuV2lpaVRlCtLiXOOtr5hSnLjEwjjF9hp1HEEkXmLJhAOAFvMbIOZBYG7gD1T1tkDfNq7fyfw7GzdP2YWMLMS734a8EHgyHyLl6XvYv8IQ6Ph+LUQstRCEFmoOTtZvWMC9wJ7AT/wgHPuqJl9Gah2zu0Bvgt838zOEGkZ3DX+fDOrBfKAoJl9BHgfcB7Y64WBH3ga+HZc35ksCTVt/QCUxqmFkJeRhplaCCILEdVRN+fcE8ATU5Z9cdL9IeBjMzy3aobNXhddibKc1bT1AVASpxaC32eU5aarhSCyABqpLEl1tq2PgM9iug7CVKvzMzUWQWQBFAiSVDVt/RTnBPHZdCeqLUxFQcbEFdhEJHqJO1FbZBo17f1xO6A8bnV+Js+eaMU5h8UxaKLlnOPZE63sWl/ELw83TSzX+ARJdWohSNKMjIW50DEQt1NOx63Oz2BoNEzXQOKvrzwaCvOfHjnE3Q9W858ffSPhry8SC7UQJGkudAwQCru4txDGZ01t7B6kMHvhl+Scr96hUX7w6nnqOge5aWMRTx9voSI/g8tW5yWsBpFYqIUgSXM2zmcYjVvtBUJTV+IOLA+PhfjG82dp7hnikzeu4wNXrKY0N53H32xkNBROWB0isVAgSNLEewzCuIr8DACauhN3YPn8xQG6Bkb5o13r2FGRT8DnY/dVFXQOjPKbk5pyRZYGBYIkTU1bH6W56XGZ1G6ykpx00vyW0MFpte39+Aw2l+VMLNtYmsNVa/J54XQb7X26rKekPgWCJE1Nez8bS7Ljvl2fz1iVl5HQwWm1FweoKMgkGHjrr9QdV6zGb8ZvT6uVIKlPgSBJc7atj42lOXOvuAAV+ZkJayEMj4Wo7xxgfdGlF/jJy0hjW3kux5p6CYU1u7ukNgWCJEVH/whdA6NsKo1/CwFgdQIHpx1p6GYs7KiaobWzoyKP/uExDp7vTEg9IgulQJCkGJ/DaNNitRAKMmnpGSKcgP/K95+L/KFfXzx9IGxblYvfZ+w92rzotYjEQoEgSXGqZZEDIT+D0ZBLyMHc6toOSnLSZ7xCW3qany1lOTx5pBldFFBSmQJBkuLR6jqy0wOLcrD14X0XJgJnsY8jhMOO6vOdVBVfevxgsu2r82joGuRoY8+i1iMSCwWCJEXtxX6qirMWba6hRF1K83RrH92Do1TN0F007vLVeRjwd78++ZbrL4ukEgWCJFxj1yCdA3P/EY1Foi6leaA2cqXY9XO0ELLTA2woyVYLQVKaAkESbvyP6Exn5cRDZtBPmt8WvYVwoLaDstx0iqKYM2l7RR6tvcO09WqQmqQmBYIk3P5zHaQHfKz2pphYDGZGfmYajYs8fUV1bSfXVxVF1fW13Zvk7mhj96LWJLJQCgRJuAO1HawvzorrRXGmU5ydPjFf0mJo6BqkoWuQ66sKo1q/ICtIRX4Gp1v7Fq0mkVgoECShOvtHONXSt6jHD8atyosEwmLNNlrtdX3tqiqK+jkbS3Oo6xhgaDS0KDWJxEKBIAn1u4OwiQiEDEZCYWrbF6eVUF3bSXbQz2XluVE/Z2NJNmNhx2sXNGpZUo8CQRLqQG0HwYCPNYWZi/5aq/IixyhOtvQuyvafOtbC6oJMflRdH/VzqkqyMeDVsxcXpSaRWCgQJKH213Zy9ZoC0vyLv+uV5qbj9xknm+MfCD1Do7T0DE07od1sMtL8VBZm8kqNAkFSjwJBEqZ/eIyjDd1cvyG6g7CxSvP7qCrOWpRAeP1CF46FdX1tLMnmUF0XgyM6jiCpRYEgCfP6hS7Gwo7r53EQNlaXlectSpdRdW0HPoO1RfPv+tpYmsNoyFF9viPudYnEQoEgCfPrY82k+Y3r1iemhQCwdVUuFzoGGBgZi+t2q2s7WZ2fSXpg/ld7W1+cRcBnvKpuI0kxCgRJiK6BER6trueKygIef6MpYa+7rTwX5+B0S/zO/R8NhXm9rpN1c0xXMZP0gJ8r1+Tzig4sS4pRIEhCPLz/AoOjId6+uTihrzt+Smg8u42ONfYwNBqOaSzFTRuLebO+m/7h+LZcRGKhQJBFNzIW5sGXa7l1Swmr8xf/dNPJ1hZlkZHmi+uB5YmxFPM8w2iyt20qZizsJrYlkgoUCLLofvFmIy09w9x9y4aEv7bfZ2xdlcupOLYQDp7vZG1RJnnejKoLcba1H7+ZTj+VlKJAkEXlnOM7vz3HlrIc3rG1NCk1bF2Vy4k4tRCccxyo7WTX+tjOlAoGfKwrzuL5k/G/QJDIQikQZFG9cvYix5p6+MytGxbtYjhzuaw8l7beYTr6R2Le1oWOAdr7htkV5YR2s9m+Oo8Tzb2cW6SpNUTmS4Egi6Z7cJT/9rPD5GUEGBoNJ+1KYVtXeQeW49BKePJIMwA3byqJeVs7KiLTYe892hzztkTiQYEgiyIcdnzukUPUdw7y8RvWJWSqiuk8vO8Cx5oiVymL9TiCc45HD9Zz3fpCNsTh4j4FWUGuXJM/ETIiyaZAkEXxtefO8MyJVr74oe0Jmdl0NrnpATLT/DEfR3i9roszrX2sK8qKW2vn/TvKOVTXRdMiX8hHJBpRBYKZ3W5mJ83sjJl9fprH083sEe/xfWZW5S0vNrPnzKzPzL425TnXmdlh7zl/b8nqYJa4e+5EK/c/dYpr1hbgT4Efq5lRnp/BkYbYrlT2aHU9aX7jisr8OFUWCQSAXx9tids2RRZqzkAwMz/wdeAOYDvwcTPbPmW1u4FO59xm4H7gPm/5EPA3wF9Os+lvAPcAW7zb7Qt5A5JajjX2cO/Dr1Gen8HuqyuTdiB5qm2rcjnc0M2FiwMLev7gSIhfvNHIzop8MtLmP13FTDaX5bC5LEfHESQlRNNCuAE445yrcc6NAD8Edk9ZZzfwoHf/x8C7zcycc/3OuReJBMMEM1sN5DnnXnHOOeAh4COxvBFJvm/+5iwf//ar+H3Gp95WRTCQOj2SV66J/Ff/+JuNC3r+3qPN9A6PLco8TLfvKGffuY64nAUlEotofmMrgbpJ39d7y6Zdxzk3BnQDs81RUOltZ7ZtyhLSPzzGQ6/UMjga4tM3V5Efw6CtxVCQFeSGqiIee72ByP8g8/PowTrWFWVRFYeDyVPdvrOcUNjx9HF1G0lyRRMI07X5p/5GRbPOgtY3s3vMrNrMqtvaNIgnVX3+p4dp6h7i49evTfj0FNH68NUVnG7tm/fB5bqOAV46c5E7r1uDbxG6wHZU5FFZkMm3X6jh4X0XknZ6rkg0gVAPrJ30/Rpgart7Yh0zCwD5wGyTtNR725ltmwA4577lnNvlnNtVWpqcka4yu2eOt/D4G428+/IytpXnJbucGX3gitUEfMZjhxrm9bz7njxB0O/jzuvWzL3yApgZf7hrLadb+2juGZr7CSKLJJpAOABsMbMNZhYE7gL2TFlnD/Bp7/6dwLNulna5c64J6DWzm7yziz4F/Hze1UvS/dNL5/jcj96gLDed25I0NUW0irKD3La1lMcPNRIOR9dt9MzxFn7xZhO3bS3hN4s4zcSn3raeoN/Hb0+pFSzJE5hrBefcmJndC+wF/MADzrmjZvZloNo5twf4LvB9MztDpGVw1/jzzawWyAOCZvYR4H3OuWPAnwPfAzKBX3k3WWKePtZCz+Aof3bbRgK+1DmIPJPdV1fw7IlWqs93csOG2ecj6hse468fO7KoYTe5e+j6qkJeqbnIe7avWpTXEpnLnIEA4Jx7AnhiyrIvTro/BHxshudWzbC8GtgZbaGSet6o6+Llsxe5YUMR65I8+Cxa77l8FZlpfh471DBnIPx/e0/S3DPEn922KSFhd8uWUl6t6eDF0+189l2bF/31RKZK/X/pJCU55/jvjx8lNyMwMbhqKchOD3DHFeU8Wl3HU8dmPqvnmeMtPPhKLZ+6aT3rYrjuwXzkZ6Zx9boCqs93cLFvOCGvKTKZAkEW5EBtJ69d6OKd28riOlBrMY2fwbNjdT7bK/L58x8cvGQeoVDY8b+eOsVnHqpm26pc/svtlyW0xlu3lDAWcnzv5dqEvq4IRNllJDLVPz5/luLs4KIM1FpsmUE/37/7Bv7NA/v57MOv8YU7LqM8P4NQ2PHjg/X89nQ7f3BtJVdWFrDn0MIGsi1UWW4GOyrz+ccXanjntlKui/G6CyLzoUCQebv/qVM8c6KV91xelrRZTGOVl5HGQ3ffyL/9p/185ZfHJ5YHfMZHr6nkunWFSZt24yNXVfCDfee556GDPPbZt7M2QV1WIraQUZvJsmvXLlddXZ3sMla8P/iHlznc0MV/ff9lZKUv7f8pQmHHDRuKMItcbvO5E61kBZP/nm7cWMRHv/4Sq/Iy+Mm/v5m8jNQa+S1Li5kddM7tmmu95O/5sqQ0dQ/yRl0XN2wsWvJhAJEQOHi+c+L7VAgDgE2lOXzzj6/jUw/sZ/fXXuIPd629ZDqQT9y4LknVyXK1NNv7kjT/9FItDsctcbhimMzu5s0l/M+PXUldxwD/++lT7D/XQXgJtehl6VEgSNSGRkP8cP8FdlTkU5gdTHY5K8JHr1nDf3z3FioKMnnsUAPff+X8gibnE4mGAkGi9utjLfQMjXF9lc58SaTinHQ+c8sG3n15GSdbeqld4DUdROaSGh2msiQ8Wl1HZUEmG0uXxqjkpWzqjKdmxq2bS3nxdDuvne+MyzWdRaZSC0Gi8g/PneHF0+1sXZW7KFNAy9yCAR9XVOZzuKGb4dFQssuRZUiBIFF5va4LB0tyINpysmt9ISOhMIdjvD60yHQUCDIn59xEN0WRDiYn1dqiLEpz0t9yqqxIvCgQZE4Haju52D/CdevUOkg2M+O69YWc7xjgbFtfssuRZUaBIHN6tLqOYMDHzsr8ZJciwDXrCvAZPFpdP/fKIvOgQJBZDY2GeOJwE1dU5BMMaHdJBbkZaWxblctPXqsnFOWV30Siod9wmdVvTrbRPxLiyrVqHaSSK9YU0NY7zPGmnmSXIsuIAkFm9Ys3GynKDrKxJCfZpcgkVcWRGVCrazuSXIksJwoEmdH3Xqpl79FmNpfm4Pdp7EEqKcgKUlmQyYFanW0k8aNAkBmdbOllNOS4Yo26i1LR9VWFHKjt0NxGEjcKBJnR4fouctIDmiYhRe2qKqK1d5gLHZrbSOJDgSDT6h8e42RLLzsr8zRVRYpq6x0G4P88cybJlchyoUCQaT1zojXSXVRZkOxSZAaluelkpvmpvdif7FJkmVAgyLR++WYjuRkB1hfrer6pymfG+uIsTYctcaNAkEv0DY/xm5Nt7KzIV3dRiqsqzqa9b5j2vuFklyLLgAJBLvHciVaGx8KaqmIJ+N14BJ1+KrFTIMglnjzSTElOUN1FS0BFQSYBn3FAA9QkDhQI8hZDoyGeOtbCxpIcdRctAQG/jzWFWRqxLHGhQJC3eP5UGyOhMDsq85JdikSpqiSLI4099A+PJbsUWeIUCPIWTx5pJjPNr7mLlpCNJTmEwo6XzrQnuxRZ4hQIMmFkLMzTx1u4fHWe5i5aQqpKsshND/DM8dZklyJLnAJBJrx0tp3eoTF2Vqi7aCkJ+Hzctq2UZ0+2Etb1ESQGCgSZ8OThZnLSA2wuU3fRUvOey8to6x3mcEN3skuRJUyBIACMhcI8dbyF37usjIBfu8VS886tZfgMnjnekuxSZAnTb74AcKC2k47+EW7fWZ7sUmQBfnWkmbVFWTx6UNdZloVTIAgAe482Ewz4eMfW0mSXIgt0eXkeTd1DNHYNJrsUWaKiCgQzu93MTprZGTP7/DSPp5vZI97j+8ysatJjX/CWnzSz909aXmtmh83skJlVx+PNyMI45/jZ6w1sKsnm54cak12OLNBl5blAZKZakYWYMxDMzA98HbgD2A583My2T1ntbqDTObcZuB+4z3vuduAuYAdwO/AP3vbGvcs5d7VzblfM70QW7HBDN92Do2yv0NxFS1lpbjpF2UEdR5AFi6aFcANwxjlX45wbAX4I7J6yzm7gQe/+j4F3m5l5y3/onBt2zp0DznjbkxSy92gzPoPLvf8wZWkyMy4vz+Xlsxc1alkWJJpAqATqJn1f7y2bdh3n3BjQDRTP8VwH/NrMDprZPTO9uJndY2bVZlbd1tYWRbkyX3uPtlBVnE1WeiDZpUiMtlfkMzIW5pdvNiW7FFmCogmE6YasTh39MtM6sz337c65a4l0RX3WzG6b7sWdc99yzu1yzu0qLdUBz3g709rHmdY+tmsw2rJQVZzFZeW5fOfFGpzTIDWZn2gCoR5YO+n7NcDUI48T65hZAMgHOmZ7rnNu/Gsr8DPUlZQUe482A7B9tQJhOTAz/vTWjZxq6eOF05rbSOYnmkA4AGwxsw1mFiRykHjPlHX2AJ/27t8JPOsi/57sAe7yzkLaAGwB9ptZtpnlAphZNvA+4Ejsb0fm69dHm7lqTT4FWcFklyJx0j8yRm5GgO/8tibZpcgSM2cgeMcE7gX2AseBHznnjprZl83sw95q3wWKzewM8Dng895zjwI/Ao4BTwKfdc6FgFXAi2b2BrAf+KVz7sn4vjWZS1P3IG/Ud/O+HRqMtpwEfD7etrGY355u52Rzb7LLkSUkqqOIzrkngCemLPvipPtDwMdmeO5Xga9OWVYDXDXfYiW+njwS6S66fWc5+2p0gZXl5IaqIl443cZ3X6zhb+/Ur5pERyOVV7AHXz5PWW66wmAZykoPcOd1a3js9UZae4eSXY4sEQqEFaqtd5jzF/vZWanBaMvVZ27ZSNg5/vbJk8kuRZYIBcIK9etjzThgh043XbZePnuRt28u4ccH69l/Tq1AmZsCYYV68kgzxdlByvMykl2KLKJ3bSujICuNv37sMKOhcLLLkRSnQFiBugZGeOXsRXZU5BOZYUSWq2DAx4eurOBUSx8PvHgu2eVIilMgrEBPH29lLOzYWanuopXg8tV5vOfyVfzvp09T1zGQ7HIkhSkQVqAnjzRRWZBJZUFmskuRBPnSh7cT8Bl/+lA1vUOjyS5HUpQCYYXpGx7jhdPtvH9HubqLVpAXTrXzsV1rOdXSy70Pv86YjifINBQIK8zTx1oYGQtzxxUanbzSbC7LYffVlTx/qo0vPX5Uk9/JJTTf8QrzzefPkp+ZxsnmXnxqIaw411cVUZaXzj8+X0NhVpDPvXerWooyQS2EFaRrYITTLX1cWZmvMFjB1hZmsWt9If/n2TP81WNHCIXVUpAItRBWkCePNBNyjivXFCS7FEkinxkfvaaS7PQAD++7QGf/CPf/0dVkpPnnfrIsa2ohrCB73mikODtIRYEGo610Zsb7d5TzgZ3l/OpIM++7/wXaeoeTXZYkmQJhhWjtGeKVmotcuaZAfcYy4ZYtpdx1/VqaugfZ/bUXOdrYneySJIkUCCvELw834RxcuUaT2clbXbmmgHtu3YQD7vzGKzz2ekOyS5IkUSCsEI+/0chl5bms0txFMo3Kwkz+zc1VlOWm8xePHOLeh1+ja2Ak2WVJgikQVoC6jgFeu9DFh66qSHYpksJyM9L4zK0bed/2VTzpHVd46liLxiusIAqEFeCtngrLAAAMbElEQVSLPz+KAej3Wubg9xnv3FbGv3vHJgD+9KFq/vi7+zje1JPkyiQRFAjL3GgoTPX5DrauyqUwO5jscmSJqCjI5P/5vS188MrVHG3s4ff//rf85x+9wYlmBcNypnEIy9xTx1roHRrjxquLkl2KLDF+n3HzphKuXlvAcyda2fNGAz95rZ5bt5Twb99exW1bSgn49T/lcqJAWOZ+8Op5CrLS2Fqem+xSZInKCgb4/SsreNdlZew/18ErNRf5k+9VU5qbzkevqeTDV1WwoyJPpzMvAwqEZexsWx8vn73I+7av0lQVErOsYIB3bivjli0lnGru5eCFLr7z2xq+9UINJTlBbtlcws2bSrhqbQGby3Lw+7TPLTUKhGXsn1+9QMBnXLe+MNmlyDIS8PnYXpHP9op8+obHONXcy+nWXp461sJjhxoBCPp9XLu+gLdtLOFtm4q5em0BwYC6l1KdAmGZGhwJ8eODdbx/Zzm5GWnJLkeWqZz0ANeuL+Ta9YWEnaO9b5iGzkHqOge50NHPvppT3P80ZKb5uXlTMe/cVso7tpaxrjgr2aXLNBQIy9QjBy7QMzTGH9+4nnPt/ckuR1YAnxlluRmU5WZwzbpIq3RwJMS59j5Ot/bxel0Xz5xoBY6ytiiTm73WwzXrClhXlKVjEClAgbAMdfaPcN+TJ9lYmk1NW59+0SRpMoP+ie4l5xwX+0Y41dpLTVs/P3+jgUeq6wAoyErjqjUF3LK5hHdsK2VLWY722yRQICxDf/fUSYbHQnzwygr9UknKMDNKctMpyU3n5k0lhJ2juXuI+s5B6jsHONbUw/On2vjqE8dZnZ/B+3eU8/tXrua6dYX4dIA6IRQIy8yxxh4e3neBGzcUU655iySF+cyoKMikoiCTGzZExsl0DYxwurWPk829/ODV83zv5VrK8zJ412WlvN07i6lIAywXjQJhGXHO8d8fP0p+ZhrvuXxVsssRmbeCrCDXVxVxfVURw6MhTjT3crihm5++1sC/7K/DDDaUZLOzIp+dlXlsXZXLxpIcKgszdZprHCgQlpEfVdex71wHX/nITo07kCUvPc3PVWsLuGptAaGwo7FrkLNtfdR3DvLCqTb2vNE4sW4w4GNdUdbEbWNpNptLc9i8KofSnHR1nUZJgbBM/OZkK1/46WE2lWYnuxSRuPP7jLVFWawt+t3pqgPDY7T2DtPeN0xb3zAX+0Y41tjDi6fbGQmFJ9bLywiwsTSHTaU5VBVnUZ6fQXl+BqvzMyjPzyQnXX8Gx+mTWAberO/i3//za6zKy+CTN65X60BWhKz0AFXpAapK3vpPkHOO3qFIWLT0DNHWN0x77zBPHWumZ2jsku3kZgSoyM9kVX4G5XnplOdlsLogk/L8DCryM1lTmEn2CgmNlfEul7GzbX38yfcOUJQd5I9vWq8LpcuKZ2bkZaaRl5nG5rKctzw2GgrTMzhKz9AY3YOj3m2E7sExzrb28fqFTvqGxi6ZKb4oO8iawkyKs4MUZgcpzAqSkx4gNyNAdnqAvIw08jIjX3MyAuSkB8gK+skOBpbUGVIKhCVszxuN/LefHibsHP/6piryNCJZZFZpfh/FOekU56TPuE4o7OgdioRF1+AoXf0jdAyM0jUwwsm+XgZGQgyMhBgZC8+4jXFmkJseIC8zjbWFWVyxJp+dlfncvKmYkllqSBYFwhI0NBriy784xsP7LnDtugLec/kqCrJ0Kp5IPPh9RkFWkIKsIOtnWS/sHCNjYYbHwgyOhhgaCTE4GmJ4LOwtDzE0GmZoNLK8rnOA/bUdhMKONL/x9s0l3LallIw0P5+4cV3C3t9sogoEM7sd+P8BP/Ad59z/O+XxdOAh4DrgIvBHzrla77EvAHcDIeA/OOf2RrNNuVR73zA/ePU8P3j1PO19I9y2pYT3bi/X6XYiSeAzIyPNT0aan/zM6FrnobCjqXuQF8+085uTbew/18HvXVbGx3atIS0Fri0xZyCYmR/4OvBeoB44YGZ7nHPHJq12N9DpnNtsZncB9wF/ZGbbgbuAHUAF8LSZbfWeM9c2V7zRUJgTTb3sO3eRfec6eP5UGyNjYd61rZTNZblsKNEZRSJLid9nrCnM4q7r13Hr5kGePNrEL95s4lRLL1/60A5u3lyS1PqiaSHcAJxxztUAmNkPgd3A5D/eu4Evefd/DHzNIif+7gZ+6JwbBs6Z2Rlve0SxzSXHOUco7Ag5x1gochsNR5qUQ6OhidvgSJiBkTEGR0MMev2RvZMOcjX3DHKhY4DGriFC4cjhraLsINeuK+CmjcWU5WoEsshSV1mYyZ+8fQPHm3p5/nQrn/jOvsg1JTYXs2t9EZetziUrzZ/Qq9JFEwiVQN2k7+uBG2daxzk3ZmbdQLG3/NUpz6307s+1zbj58Nde5GxrH/C768w7747DTdyPhpu0kfHnhp0jHIcL2KcHfGSm+cnNCFCUHWRzaQ5luRlUlWRH3SQVkaXDzNhekcdff/Byvv1CDT871MDfPnnyLesEfEZ6wEf1X7+XzODinkUYTSBM10E99c/fTOvMtHy6yJv2T6qZ3QPc433bZ2Ynp1svQUqA9iS+/kxStS5QbQuVqrWlal2wzGvL+h8xvf5sx8cnRBMI9cDaSd+vARpnWKfezAJAPtAxx3Pn2iYAzrlvAd+Kos5FZ2bVzrldya5jqlStC1TbQqVqbalaF6i2eIimc+oAsMXMNphZkMhB4j1T1tkDfNq7fyfwrHPOecvvMrN0M9sAbAH2R7lNERFJoDlbCN4xgXuBvUROEX3AOXfUzL4MVDvn9gDfBb7vHTTuIPIHHm+9HxE5WDwGfNY5FwKYbpvxf3siIhKtqMYhOOeeAJ6YsuyLk+4PAR+b4blfBb4azTaXgJTouppGqtYFqm2hUrW2VK0LVFvMzM3nFBsREVm2kj80TkREUoICYRpmts3MDk269ZjZX5jZl8ysYdLyDySongfMrNXMjkxaVmRmT5nZae9robfczOzvzeyMmb1pZtcmobb/aWYnvNf/mZkVeMurzGxw0uf3zQTXNePPz8y+4H1mJ83s/YtV1yy1PTKprlozO+QtT9hn5r3eWjN7zsyOm9lRM/uP3vKk7m+z1JUK+9pMtaXE/jYvzjndZrkROejdTOQ83i8Bf5mEGm4DrgWOTFr2t8DnvfufB+7z7n8A+BWRMSA3AfuSUNv7gIB3/75JtVVNXi8JdU378wO2A28A6cAG4CzgT2RtUx7/O+CLif7MvNdbDVzr3c8FTnmfT1L3t1nqSoV9babaUmJ/m89NLYS5vRs465w7n6wCnHMvEDl7a7LdwIPe/QeBj0xa/pCLeBUoMLPViazNOfdr59z4lUheJTLOJKFm+MxmMjHFinPuHDB5ipWE1mZmBvwh8C+L9fqzcc41Oede8+73AseJzC6Q1P1tprpSZF+b6TObSUL3t/lQIMztLt76y3mv1zx9YLzZnCSrnHNNENkhgTJv+XRTjcy2cy62PyHyH+S4DWb2upk9b2a3JqGe6X5+qfSZ3Qq0OOdOT1qWlM/MzKqAa4B9pND+NqWuyZK+r01TW6rvb2+hQJiFRQbNfRh41Fv0DWATcDXQRKRpn2qimWokIczsr4iMP/lnb1ETsM45dw3wOeBhM8tLYEkz/fxS5jMDPs5b/wFJymdmZjnAT4C/cM71zLbqNMsW7bObqa5U2NemqW0p7G9voUCY3R3Aa865FgDnXItzLuScCwPfJrnNvJbxprn3tdVbHs1UI4vOzD4NfBD4pPM6Tr0m8kXv/kEifadbZ95KfM3y80uVzywA/AHwyPiyZHxmZpZG5A/bPzvnfuotTvr+NkNdKbGvTVdbqu9v01EgzO4t/61N6Rv9KHDkkmckzuTpQj4N/HzS8k95Z3/cBHSPN/UTxSIXP/qvwIedcwOTlpda5PoamNlGIlOZ1CSwrpl+fjNNsZJo7wFOOOfqxxck+jPzjmF8FzjunPtfkx5K6v42U12psK/NUluq72+XSvZR7VS9AVlErv6WP2nZ94HDwJtEfqirE1TLvxBpco4S+e/ibiLTiz8DnPa+FnnrGpGLD531at2VhNrOEOkjPeTdvumt+6+Ao0TOsHgN+FCC65rx5wf8lfeZnQTuSPRn5i3/HvDvpqybsM/Me71biHRfvDnp5/eBZO9vs9SVCvvaTLWlxP42n5tGKouICKAuIxER8SgQREQEUCCIiIhHgSAiIoACQUREPAoEEREBFAgiIuJRIIiICAD/F0gOuLOA5UcaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.distplot(S_.flatten(),  bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 566,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv2.imshow('orig', S/255)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()\n",
    "\n",
    "cv2.imshow('GC', S_/255)\n",
    "cv2.waitKey(0)\n",
    "cv2.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}