demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (c) Nanjing University, Vision Lab.\n",
    "# Last update: 2020.12.09"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import time\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import matplotlib.pylab  as plt\n",
    "import pandas as pd\n",
    "tf.enable_eager_execution()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.13.1\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "print(tf.__version__)\n",
    "print(tf.test.is_gpu_available())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### set gpu/cpu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# os.environ['CUDA_VISIBLE_DEVICES']=\"\" # only cpu\n",
    "\n",
    "config = tf.ConfigProto()\n",
    "config.gpu_options.per_process_gpu_memory_fraction = 0.5\n",
    "config.gpu_options.allow_growth = True\n",
    "config.log_device_placement=True\n",
    "# config.device_count={'gpu':0}\n",
    "sess = tf.Session(config=config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Compression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "###  Pre-process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from process import preprocess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file = 'testdata/8iVFB/longdress_vox10_1300.ply'\n",
    "scale = 1.\n",
    "cube_size = 64\n",
    "min_num = 64\n",
    "res = 1024"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Preprocess =====\n",
      "Scaling: 0.0s\n",
      "Partition: 10.1974s\n",
      "Voxelization: 0.1488s\n",
      "cubes shape: (202, 64, 64, 64, 1)\n",
      "points numbers (sum/mean/max/min): 857600 4246 11450 78\n"
     ]
    }
   ],
   "source": [
    "cubes, cube_positions, points_numbers = preprocess(input_file, scale, cube_size, min_num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "point_numbers (min/max/mean):  78 11450 4245\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAACcCAYAAACwV5WjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPZklEQVR4nO3dfawldX3H8fengGgBYWG3FOXhgqKGGl3xipoaYqXyoEa0sXaJqcQntNXERmOzaFs1sRVqWh8i1VABMVWgYqmU1fKgNJio6F0FXEXcha5lEdjFlSdrVfDbP+a3y+HmPuyce3bvw75fycmZ+c2cmd/5yqyfOzPnN6kqJEmStON+a747IEmStNgYoCRJknoyQEmSJPVkgJIkSerJACVJktSTAUqSJKknA5SkkUjyYJKj5rsfo5Skkjx5vvshaeExQEkaiarat6pu25F1DSaSFjsDlCTtZEn2nO8+SBotA5Sk7ZJsTHJmkh8k+VmSC5I8dmD5m5JsSLI1yeVJnjCwbPtZpSSfTnJOkjVJHkhyfZIntWXXtY/c2C77/UmS5UmuSHJv2/bXkkz571Pbz1uSrG/rn5Mkbdn7kvzLwLpjbf092/x/JflAkq+3ff9HkoOSfDbJ/Um+nWRs0i5fkuS2JPck+dBgv5K8PsnNrVZXJjliUj/fmmQ9sD6dDyfZ3Pb1vSRPH+5/KUnzzQAlabLXACcBTwKeAvwVQJIXAR8EXg0cAvwYuHiG7awC3g8sAzYAfwtQVce35c9sl/0uAd4JbAJWAAcD7wZmes7Uy4DnAM9o/Tmpx/dbBfwp8MT2Hb8BXAAcCNwMvHfS+q8ExoFjgVOB1wMkObX1849av78GXDTps68AngscA5wIHE9X0/1bv3/ao9+SFhADlKTJPl5Vt1fVVrrQc1prfw1wflV9p6p+CZwJPH+KMzbbXFZV36qqh4DPAitn2Oev6ULZEVX166r6Ws38oM6zqureqvof4NpZtj3ZBVV1a1XdB3wZuLWqrmn9/DzwrEnrn11VW9u+PsIj9XgL8MGqurl99u+AlYNnodryrVX1i/Yd9wOeBqR97s4e/Za0gBigJE12+8D0j4Ftl+me0OYBqKoH6c6gPHGa7dw1MP2/wL4z7PNDdGeprmqXy1bP0sc+257s7oHpX0wxP3lb09XjCOCj7TLivcBWIDy6Hts/W1VfBT4OnANsTnJuksf36LekBcQAJWmywwamDwd+0qZ/QhcaAEiyD3AQcMdcd1hVD1TVO6vqKODlwDuSnDDEpn4O/PbA/O/OtW9MX4/bgTdX1QEDr8dV1dcH1n/UWbSq+lhVPZvukt5TgHeNoH+S5oEBStJkb01yaJIDgfcAl7T2i4DXJVmZZG+6S1bXV9XGIfZxN7B9zKgkL0vy5HYz+H3Aw8BvhtjuDcDxSQ5Psj/dZca5eleSZUkOA97OI/X4JHBmkt8DSLJ/kj+ebiNJnpPkuUn2ogt6/8dw31HSAmCAkjTZ54CrgNuAW4EPAFTVNcBfA18A7qS7AXvVkPt4H3Bhu/z1auBo4BrgQbqbuv+pqq7tu9Gqupou4NwErAWuGLJ/g77YtnUDsAY4r+3rMuBs4OIk9wPrgFNm2M7jgX8GfkZ3KfCndJcuJS1Cmfk+TUm7kyQbgTe2sCRJmoZnoCRJknoyQEmSJPXkJTxJkqSePAMlSZLUkwFKkiSpp136hPDly5fX2NjYrtylJEnSUNauXXtPVa2YatmsASrJ+XQP7txcVU9vbe8D3gRsaau9u6q+NNu2xsbGmJiY2NF+S5IkzZskP55u2Y5cwvs0cPIU7R+uqpXtNWt4kiRJWipmDVBVdR3dQzIlSZLE3G4if1uSm5Kcn2TZyHokSZK0wA0boD5B9xyslXTPxPqH6VZMckaSiSQTW7ZsmW61kRtbvYax1Wt22f6GtVj6KUmSHjFUgKqqu6vq4ar6Dd3DMY+bYd1zq2q8qsZXrJjyRnZJkqRFZagAleSQgdlX0j2FXJIkabewI8MYXAS8EFieZBPwXuCFSVYCBWwE3rzzuihJkrSwzBqgquq0KZrP2wl9kSRJWhR8lIskSVJPBihJkqSeDFCSJEk9GaAkSZJ6MkBJkiT1ZICSJEnqyQAlSZLUkwFKkiSpJwOUJElSTwYoSZKkngxQkiRJPRmgJEmSejJASZIk9WSAkiRJ6skAJUmS1JMBSpIkqScD1C4ytnoNY6vXzHc3JEnSCBigJEmSejJASZIk9WSAkiRJ6skAJUmS1JMBSpIkqScDlCRJUk8GqF3M4QwkSVr8DFCSJEk9GaAkSZJ6mjVAJTk/yeYk6wbaDkxydZL17X3Zzu2mJEnSwrEjZ6A+DZw8qW018JWqOhr4SpuXJEnaLcwaoKrqOmDrpOZTgQvb9IXAK0bbLUmSpIVr2HugDq6qO9v0XcDBI+qPJEnSgjfnm8irqoCabnmSM5JMJJnYsmXLXHe35DnMgSRJC9+wAeruJIcAtPfN061YVedW1XhVja9YsWLI3UmSJC0cwwaoy4HT2/TpwBdH0x1JkqSFb0eGMbgI+Abw1CSbkrwBOAt4cZL1wB+2eUmSpN3CnrOtUFWnTbPohBH3RZIkaVFwJHJJkqSeDFCSJEk9zXoJb7GZbgiAbe0bz3rpUMtn21/fz0mSpMXLM1CSJEk9GaAkSZJ6MkBJkiT1ZICSJEnqyQAlSZLUkwFKkiSpJwPUAjW2es20QzJIkqT5ZYCSJEnqyQAlSZLUkwFKkiSpJwOUJElSTwYoSZKkngxQk/jrN0mSNBsDlCRJUk8GKEmSpJ4MUJIkST0ZoCRJknoyQEmSJPVkgJIkSeppz/nuwEK3bUiDjWe9dE6flyRJS4dnoCRJknoyQEmSJPU0p0t4STYCDwAPAw9V1fgoOiVJkrSQjeIeqD+oqntGsB1JkqRFwUt4kiRJPc01QBVwVZK1Sc4YRYckSZIWurlewntBVd2R5HeAq5P8sKquG1yhBaszAA4//PA57m7+zHU4g+m2J0mSFp85nYGqqjva+2bgMuC4KdY5t6rGq2p8xYoVc9mdJEnSgjB0gEqyT5L9tk0DJwLrRtUxSZKkhWoul/AOBi5Lsm07n6uq/xxJryRJkhawoQNUVd0GPHOEfZEkSVoUHMZAkiSpJwOUJElST6MYiXxRmTx8wKiHJRjV9mbb/s7enyRJmp5noCRJknoyQEmSJPVkgJIkSerJACVJktSTAUqSJKknA5QkSVJPu90wBpNNHtZg2PZh1xt2/cmfmzycgcMcSJK083gGSpIkqScDlCRJUk8GKEmSpJ4MUJIkST0ZoCRJknoyQEmSJPW05IcxGHZ4gIW6v1ENr+DwBpIkDc8zUJIkST0ZoCRJknoyQEmSJPVkgJIkSerJACVJktSTAUqSJKmnVNUu29n4+HhNTEzs1H3s6mELFrttwxn0qZtDIEiSdgdJ1lbV+FTLPAMlSZLUkwFKkiSppzkFqCQnJ7klyYYkq0fVKUmSpIVs6ACVZA/gHOAU4BjgtCTHjKpjkiRJC9VczkAdB2yoqtuq6lfAxcCpo+mWJEnSwjWXAPVE4PaB+U2tTZIkaUnbc2fvIMkZwBlt9sEkt+zkXS4H7tnJ+1gycvb2yR2u28Bndnf+tzYc69afNevPmg3Huj3aEdMtmEuAugM4bGD+0Nb2KFV1LnDuHPbTS5KJ6cZs0PSsW3/WbDjWrT9r1p81G45123FzuYT3beDoJEcmeQywCrh8NN2SJElauIY+A1VVDyV5G3AlsAdwflV9f2Q9kyRJWqDmdA9UVX0J+NKI+jIqu+xy4RJj3fqzZsOxbv1Zs/6s2XCs2w7apc/CkyRJWgp8lIskSVJPSypA7e6PlklyfpLNSdYNtB2Y5Ook69v7staeJB9rtbopybEDnzm9rb8+yekD7c9O8r32mY8lya79hqOX5LAk1yb5QZLvJ3l7a7duM0jy2CTfSnJjq9v7W/uRSa5v3/WS9gMTkuzd5je05WMD2zqztd+S5KSB9iV5PCfZI8l3k1zR5q3ZDJJsbMfPDUkmWpvH5yySHJDk0iQ/THJzkudbtxGrqiXxoruR/VbgKOAxwI3AMfPdr11cg+OBY4F1A21/D6xu06uBs9v0S4AvAwGeB1zf2g8Ebmvvy9r0srbsW23dtM+eMt/feQQ1OwQ4tk3vB/yI7tFE1m3mugXYt03vBVzfvuO/Aqta+yeBP2vTfw58sk2vAi5p08e0Y3Vv4Mh2DO+xlI9n4B3A54Ar2rw1m7leG4Hlk9o8Pmev24XAG9v0Y4ADrNtoX0vpDNRu/2iZqroO2Dqp+VS6A4n2/oqB9s9U55vAAUkOAU4Crq6qrVX1M+Bq4OS27PFV9c3qjp7PDGxr0aqqO6vqO236AeBmuhH1rdsM2vd/sM3u1V4FvAi4tLVPrtu2el4KnND+Yj0VuLiqfllV/w1soDuWl+TxnORQ4KXAp9p8sGbD8PicQZL96f6gPg+gqn5VVfdi3UZqKQUoHy0ztYOr6s42fRdwcJuerl4ztW+aon3JaJdInkV3NsW6zaJdiroB2Ez3D+utwL1V9VBbZfC7bq9PW34fcBD967nYfQT4S+A3bf4grNlsCrgqydp0T7YAj8/ZHAlsAS5ol4s/lWQfrNtILaUApVm0vxT82eUUkuwLfAH4i6q6f3CZdZtaVT1cVSvpnkJwHPC0+e3RwpbkZcDmqlo7331ZZF5QVccCpwBvTXL84EKPzyntSXc7xyeq6lnAz+ku2W1n3eZuKQWoHXq0zG7o7na6lfa+ubVPV6+Z2g+don3RS7IXXXj6bFX9W2u2bjuoXRq4Fng+3an/bePLDX7X7fVpy/cHfkr/ei5mvw+8PMlGustrLwI+ijWbUVXd0d43A5fRhXWPz5ltAjZV1fVt/lK6QGXdRmgpBSgfLTO1y4Ftv5w4HfjiQPtr268vngfc107tXgmcmGRZ+4XGicCVbdn9SZ7X7sN47cC2Fq32Xc4Dbq6qfxxYZN1mkGRFkgPa9OOAF9PdP3Yt8Kq22uS6bavnq4Cvtr+ALwdWpfvF2ZHA0XQ3py6547mqzqyqQ6tqjO77fLWqXoM1m1aSfZLst22a7rhah8fnjKrqLuD2JE9tTScAP8C6jdbOujt9Pl50vyT4Ed29GO+Z7/7Mw/e/CLgT+DXdXyBvoLtn4ivAeuAa4MC2boBzWq2+B4wPbOf1dDembgBeN9A+TveP163Ax2kDsS7mF/ACutPYNwE3tNdLrNusdXsG8N1Wt3XA37T2o+j+z3wD8Hlg79b+2Da/oS0/amBb72m1uYWBX/Is5eMZeCGP/ArPmk1fp6Pofk14I/D9bd/J43OHarcSmGjH6L/T/YrOuo3w5UjkkiRJPS2lS3iSJEm7hAFKkiSpJwOUJElSTwYoSZKkngxQkiRJPRmgJEmSejJASZIk9WSAkiRJ6un/AWLiJzdVczBBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(10,2))\n",
    "print('point_numbers (min/max/mean): ', points_numbers.min(), points_numbers.max(), \n",
    "      int(points_numbers.mean()))\n",
    "plt.title('points numbers')\n",
    "n, bins, patches = plt.hist(points_numbers, bins=256, range=(1, pow(2, 16)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from transform import compress_factorized, decompress_factorized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import models.model_voxception as model\n",
    "# import models.model_simple as model\n",
    "# import models.model_quach as model\n",
    "# import models.model_baseline as model\n",
    "ckpt_dir = './checkpoints/factorized/a2b3/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Compress =====\n",
      "WARNING:tensorflow:From /home/ubuntu/anaconda3/envs/tf13/lib/python3.6/site-packages/tensorflow/python/ops/resource_variable_ops.py:642: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n",
      "Analysis Transform: 7.4024s\n",
      "Entropy Encode: 0.1356s\n",
      "7.725878000259399\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "strings, min_v, max_v, shape = compress_factorized(cubes, model, ckpt_dir)\n",
    "print(time.time() - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2 2 [202  16  16  16  16]\n"
     ]
    }
   ],
   "source": [
    "print(min_v.numpy(), max_v.numpy(), shape.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./compressed/ longdress_vox10_1300\n",
      "===== Write binary files =====\n",
      "Total file size (Bytes): 13960\n",
      "Strings (Bytes): 13466\n",
      "Numbers of points (Bytes): 404\n",
      "Positions of cubes (Bytes): 90\n"
     ]
    }
   ],
   "source": [
    "from dataprocess.inout_bitstream import write_binary_files_factorized\n",
    "filename = os.path.split(input_file)[-1][:-4]\n",
    "rootdir = './compressed/'\n",
    "print(rootdir, filename)\n",
    "\n",
    "bytes_strings, bytes_pointnums, bytes_cubepos = write_binary_files_factorized(\n",
    "    filename, strings.numpy(), points_numbers, cube_positions,\n",
    "    min_v.numpy(), max_v.numpy(), shape.numpy(), rootdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'longdress_vox10_1300.cubepos\\n'\n",
      "b'longdress_vox10_1300_cubepos.ply\\n'\n",
      "b'longdress_vox10_1300.pointnums\\n'\n",
      "b'longdress_vox10_1300.strings\\n'\n",
      "b'longdress_vox10_1300.strings_head\\n'\n",
      "b'longdress_vox10_1300.strings_hyper\\n'\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "subp=subprocess.Popen(\"ls ./compressed/\", \n",
    "                      shell=True, stdout=subprocess.PIPE)\n",
    "c=subp.stdout.readline()\n",
    "while c:\n",
    "    print (c)\n",
    "    c=subp.stdout.readline()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Decompress"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Read binary files =====\n"
     ]
    }
   ],
   "source": [
    "from dataprocess.inout_bitstream import read_binary_files_factorized\n",
    "strings_d, points_numbers_d, cube_positions_d, min_v_d, max_v_d, shape_d = \\\n",
    "                                            read_binary_files_factorized(filename, rootdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### decoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Decompress =====\n",
      "Entropy Decode: 0.1505s\n",
      "Synthesis Transform: 6.7069s\n"
     ]
    }
   ],
   "source": [
    "from transform import decompress_factorized\n",
    "cubes_d = decompress_factorized(strings_d, min_v_d, max_v_d, shape_d, model, ckpt_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# _ = plt.hist(cubes_d.numpy().reshape(-1), bins=512)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### post-process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Post process =====\n",
      "Classify and extract points: 0.649s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.2594s\n"
     ]
    }
   ],
   "source": [
    "from process import postprocess\n",
    "rho = 1.0\n",
    "output_file = filename + \"_rec.ply\"\n",
    "postprocess(output_file, cubes_d.numpy(), \n",
    "            points_numbers_d, cube_positions_d, scale, cube_size, rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### bpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_points_number(filedir):\n",
    "    plyfile = open(filedir)\n",
    "\n",
    "    line = plyfile.readline()\n",
    "    while line.find(\"element vertex\") == -1:\n",
    "        line = plyfile.readline()\n",
    "    number = int(line.split(' ')[-1][:-1])\n",
    "    \n",
    "    return number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "857966\n"
     ]
    }
   ],
   "source": [
    "N = get_points_number(input_file)\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bpp:  0.1302\n",
      "strings:  0.1256 points numbers:  0.0038 cubes positions:  0.0008\n"
     ]
    }
   ],
   "source": [
    "print(\"bpp: \", round(8*(bytes_strings + bytes_pointnums + bytes_cubepos)/float(N), 4))\n",
    "print(\"strings: \", round(8*bytes_strings/float(N),4), \n",
    "      \"points numbers: \", round(8*bytes_pointnums/float(N),4), \n",
    "      \"cubes positions: \", round(8*bytes_cubepos/float(N),4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### pc_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'PCC quality measurement software, version 0.13.4\\n'\n",
      "b'\\n'\n",
      "b'infile1:        testdata/8iVFB/longdress_vox10_1300.ply\\n'\n",
      "b'infile2:        longdress_vox10_1300_rec.ply\\n'\n",
      "b'normal1:        testdata/8iVFB/longdress_vox10_1300.ply\\n'\n",
      "b'singlePass:     0\\n'\n",
      "b'hausdorff:      1\\n'\n",
      "b'color:          0\\n'\n",
      "b'lidar:          0\\n'\n",
      "b'resolution:     1023\\n'\n",
      "b'dropDuplicates: 2\\n'\n",
      "b'neighborsProc:  1\\n'\n",
      "b'averageNormals: 1\\n'\n",
      "b'nbThreads:      1\\n'\n",
      "b'\\n'\n",
      "b'Verifying if the data is loaded correctly.. The last point is: 256 902 320\\n'\n",
      "b'Reading file 1 done.\\n'\n",
      "b'Normals loaded from file1.\\n'\n",
      "b'Verifying if the data is loaded correctly.. The last point is: 223 960 320\\n'\n",
      "b'Reading file 2 done.\\n'\n",
      "b'Imported intrinsic resoluiton: 1023\\n'\n",
      "b'Peak distance for PSNR: 1023\\n'\n",
      "b'Point cloud sizes for org version, dec version, and the scaling ratio: 857966, 857600, 0.999573\\n'\n",
      "b'Normals prepared.\\n'\n",
      "b'\\n'\n",
      "b'1. Use infile1 (A) as reference, loop over A, use normals on B. (A->B).\\n'\n",
      "b'   mse1      (p2point): 0.45654\\n'\n",
      "b'   mse1,PSNR (p2point): 68.3739\\n'\n",
      "b'   mse1      (p2plane): 0.144076\\n'\n",
      "b'   mse1,PSNR (p2plane): 73.3828\\n'\n",
      "b'   h.       1(p2point): 107\\n'\n",
      "b'   h.,PSNR  1(p2point): 44.6749\\n'\n",
      "b'   h.       1(p2plane): 15.7771\\n'\n",
      "b'   h.,PSNR  1(p2plane): 52.9885\\n'\n",
      "b'2. Use infile2 (B) as reference, loop over B, use normals on A. (B->A).\\n'\n",
      "b'   mse2      (p2point): 0.311119\\n'\n",
      "b'   mse2,PSNR (p2point): 70.0395\\n'\n",
      "b'   mse2      (p2plane): 0.178194\\n'\n",
      "b'   mse2,PSNR (p2plane): 72.4598\\n'\n",
      "b'   h.       2(p2point): 98\\n'\n",
      "b'   h.,PSNR  2(p2point): 45.0565\\n'\n",
      "b'   h.       2(p2plane): 97.345\\n'\n",
      "b'   h.,PSNR  2(p2plane): 45.0856\\n'\n",
      "b'3. Final (symmetric).\\n'\n",
      "b'   mseF      (p2point): 0.45654\\n'\n",
      "b'   mseF,PSNR (p2point): 68.3739\\n'\n",
      "b'   mseF      (p2plane): 0.178194\\n'\n",
      "b'   mseF,PSNR (p2plane): 72.4598\\n'\n",
      "b'   h.        (p2point): 107\\n'\n",
      "b'   h.,PSNR   (p2point): 44.6749\\n'\n",
      "b'   h.        (p2plane): 97.345\\n'\n",
      "b'   h.,PSNR   (p2plane): 45.0856\\n'\n",
      "b'Job done! 11.576 seconds elapsed (excluding the time to load the point clouds).\\n'\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "\n",
    "subp=subprocess.Popen(\"myutils/pc_error_d\" \n",
    "                      + ' -a ' + input_file\n",
    "                      + ' -b ' + output_file\n",
    "                      +' -n ' + input_file \n",
    "                      + ' --hausdorff=1 '\n",
    "                      + \" -r \"+str(res-1), \n",
    "                      shell=True, stdout=subprocess.PIPE)\n",
    "\n",
    "c=subp.stdout.readline()\n",
    "while c:\n",
    "    print(c)\n",
    "    c=subp.stdout.readline()\n",
    "\n",
    "print(subp.returncode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import open3d as o3d\n",
    "pcd = o3d.io.read_point_cloud(input_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcd_rec = o3d.io.read_point_cloud(output_file)\n",
    "o3d.visualization.draw_geometries([pcd_rec])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Model with Hyperprior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_file = 'testdata/8iVFB/longdress_vox10_1300.ply'\n",
    "scale = 1.\n",
    "cube_size = 64\n",
    "min_num = 64\n",
    "res = 1024"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Preprocess =====\n",
      "Scaling: 0.0s\n",
      "Partition: 10.4049s\n",
      "Voxelization: 0.4471s\n",
      "cubes shape: (202, 64, 64, 64, 1)\n",
      "points numbers (sum/mean/max/min): 857600 4246 11450 78\n"
     ]
    }
   ],
   "source": [
    "from process import preprocess\n",
    "cubes, cube_positions, points_numbers = preprocess(input_file, scale, cube_size, min_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## compress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "import models.model_voxception as model\n",
    "from transform import compress_hyper\n",
    "ckpt_dir = './checkpoints/hyper/a0.75b3/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Compress =====\n",
      "Analysis Transform: 5.839s\n",
      "Hyper Encoder: 0.18s\n",
      "Quantize hyperprior.\n",
      "Hyper Decoder: 0.3564s\n",
      "Entropy Encode (Hyper): 0.0091s\n",
      "Entropy Encode: 1.3656s\n",
      "Entropy Decode: 1.3281s\n",
      "Synthesis Transform: 6.6119s\n",
      "15.9921\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "y_strings, y_min_vs, y_max_vs, y_shape, z_strings, z_min_v, z_max_v, z_shape, x_ds = \\\n",
    "compress_hyper(cubes, model, ckpt_dir, True)\n",
    "print(round(time.time() - start, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Write binary files =====\n",
      "Total file size (Bytes): 12148\n",
      "Strings (Bytes): 7128\n",
      "Strings head (Bytes): 416\n",
      "Strings hyper (Bytes): 4110\n",
      "Numbers of points (Bytes): 404\n",
      "Positions of cubes (Bytes): 90\n"
     ]
    }
   ],
   "source": [
    "from dataprocess.inout_bitstream import write_binary_files_hyper\n",
    "filename = os.path.split(input_file)[-1][:-4]\n",
    "rootdir = './compressed/'\n",
    "\n",
    "bytes_strings, bytes_strings_head, bytes_strings_hyper, bytes_pointnums, bytes_cubepos = \\\n",
    "    write_binary_files_hyper(\n",
    "    filename, y_strings.numpy(), z_strings.numpy(), points_numbers, cube_positions,\n",
    "    y_min_vs.numpy(), y_max_vs.numpy(), y_shape.numpy(), \n",
    "    z_min_v.numpy(), z_max_v.numpy(), z_shape.numpy(), rootdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## decompress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Read binary files =====\n"
     ]
    }
   ],
   "source": [
    "from dataprocess.inout_bitstream import read_binary_files_hyper\n",
    "y_strings_d, z_strings_d, points_numbers_d, cube_positions_d, \\\n",
    "y_min_vs_d, y_max_vs_d, y_shape_d, z_min_v_d, z_max_v_d, z_shape_d = \\\n",
    "read_binary_files_hyper(filename, rootdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== Decompress =====\n",
      "Entropy Decoder (Hyper): 0.026s\n",
      "Hyper Decoder: 0.3464s\n",
      "Entropy Decoder: 1.3539s\n",
      "Synthesis Transform: 6.6596s\n"
     ]
    }
   ],
   "source": [
    "from transform import decompress_hyper\n",
    "cubes_d = decompress_hyper(y_strings_d, y_min_vs_d.astype('int32'), y_max_vs_d.astype('int32'), y_shape_d, \n",
    "                     z_strings_d, z_min_v_d, z_max_v_d, z_shape_d, model, ckpt_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error on gpu: 1.1444092e-05\n"
     ]
    }
   ],
   "source": [
    "# check\n",
    "print(\"error on gpu:\", tf.reduce_max(tf.abs(cubes_d - x_ds)).numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "longdress_vox10_1300_rec.ply\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.359s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.3553s\n"
     ]
    }
   ],
   "source": [
    "from process import postprocess\n",
    "# cubes_d = x_ds\n",
    "rho = 1.0\n",
    "output_file = filename + \"_rec.ply\"\n",
    "print(output_file)\n",
    "\n",
    "postprocess(output_file, cubes_d.numpy(), \n",
    "            points_numbers_d, cube_positions_d, scale, cube_size, rho)\n",
    "# cheat\n",
    "# postprocess(output_file, x_ds, points_numbers_d, cube_positions_d, scale, cube_size, rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "857966.0\n",
      "bpp:  0.1133\n",
      "strings:  0.0665 strings_hyper:  0.0383 \n",
      " strings_head:  0.0039 points numbers:  0.0038 cubes positions:  0.0008\n"
     ]
    }
   ],
   "source": [
    "def get_points_number(filedir):\n",
    "    plyfile = open(filedir)\n",
    "\n",
    "    line = plyfile.readline()\n",
    "    while line.find(\"element vertex\") == -1:\n",
    "        line = plyfile.readline()\n",
    "    number = float(line.split(' ')[-1][:-1])\n",
    "    \n",
    "    return number\n",
    "\n",
    "N = get_points_number(input_file)\n",
    "print(N)\n",
    "\n",
    "print(\"bpp: \", \n",
    "      round(8*(bytes_strings \n",
    "               + bytes_strings_head + \n",
    "               bytes_strings_hyper + \n",
    "               bytes_pointnums + \n",
    "               bytes_cubepos)/N,4))\n",
    "print(\"strings: \", round(8*bytes_strings/N,4), \n",
    "      \"strings_hyper: \", round(8*bytes_strings_hyper/N,4), '\\n', \n",
    "      \"strings_head: \", round(8*bytes_strings_head/N,4), \n",
    "      \"points numbers: \", round(8*bytes_pointnums/N,4), \n",
    "      \"cubes positions: \", round(8*bytes_cubepos/N,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'PCC quality measurement software, version 0.13.4\\n'\n",
      "b'\\n'\n",
      "b'infile1:        testdata/8iVFB/longdress_vox10_1300.ply\\n'\n",
      "b'infile2:        longdress_vox10_1300_rec.ply\\n'\n",
      "b'normal1:        testdata/8iVFB/longdress_vox10_1300.ply\\n'\n",
      "b'singlePass:     0\\n'\n",
      "b'hausdorff:      1\\n'\n",
      "b'color:          0\\n'\n",
      "b'lidar:          0\\n'\n",
      "b'resolution:     1023\\n'\n",
      "b'dropDuplicates: 2\\n'\n",
      "b'neighborsProc:  1\\n'\n",
      "b'averageNormals: 1\\n'\n",
      "b'nbThreads:      1\\n'\n",
      "b'\\n'\n",
      "b'Verifying if the data is loaded correctly.. The last point is: 256 902 320\\n'\n",
      "b'Reading file 1 done.\\n'\n",
      "b'Normals loaded from file1.\\n'\n",
      "b'Verifying if the data is loaded correctly.. The last point is: 220 960 323\\n'\n",
      "b'Reading file 2 done.\\n'\n",
      "b'Imported intrinsic resoluiton: 1023\\n'\n",
      "b'Peak distance for PSNR: 1023\\n'\n",
      "b'Point cloud sizes for org version, dec version, and the scaling ratio: 857966, 857600, 0.999573\\n'\n",
      "b'Normals prepared.\\n'\n",
      "b'\\n'\n",
      "b'1. Use infile1 (A) as reference, loop over A, use normals on B. (A->B).\\n'\n",
      "b'   mse1      (p2point): 0.531367\\n'\n",
      "b'   mse1,PSNR (p2point): 67.7148\\n'\n",
      "b'   mse1      (p2plane): 0.190819\\n'\n",
      "b'   mse1,PSNR (p2plane): 72.1625\\n'\n",
      "b'   h.       1(p2point): 94\\n'\n",
      "b'   h.,PSNR  1(p2point): 45.2374\\n'\n",
      "b'   h.       1(p2plane): 22.7506\\n'\n",
      "b'   h.,PSNR  1(p2plane): 51.3988\\n'\n",
      "b'2. Use infile2 (B) as reference, loop over B, use normals on A. (B->A).\\n'\n",
      "b'   mse2      (p2point): 0.388849\\n'\n",
      "b'   mse2,PSNR (p2point): 69.0709\\n'\n",
      "b'   mse2      (p2plane): 0.2414\\n'\n",
      "b'   mse2,PSNR (p2plane): 71.1413\\n'\n",
      "b'   h.       2(p2point): 37\\n'\n",
      "b'   h.,PSNR  2(p2point): 49.2867\\n'\n",
      "b'   h.       2(p2plane): 35.2584\\n'\n",
      "b'   h.,PSNR  2(p2plane): 49.4961\\n'\n",
      "b'3. Final (symmetric).\\n'\n",
      "b'   mseF      (p2point): 0.531367\\n'\n",
      "b'   mseF,PSNR (p2point): 67.7148\\n'\n",
      "b'   mseF      (p2plane): 0.2414\\n'\n",
      "b'   mseF,PSNR (p2plane): 71.1413\\n'\n",
      "b'   h.        (p2point): 94\\n'\n",
      "b'   h.,PSNR   (p2point): 45.2374\\n'\n",
      "b'   h.        (p2plane): 35.2584\\n'\n",
      "b'   h.,PSNR   (p2plane): 49.4961\\n'\n",
      "b'Job done! 11.667 seconds elapsed (excluding the time to load the point clouds).\\n'\n",
      "None\n"
     ]
    }
   ],
   "source": [
    "import subprocess\n",
    "\n",
    "input_file_n = input_file[:-4]+'.ply'\n",
    "subp=subprocess.Popen(\"myutils/pc_error_d\" \n",
    "                      + ' -a ' + input_file\n",
    "                      + ' -b ' + output_file\n",
    "                      +' -n ' + input_file \n",
    "                      + ' --hausdorff=1 '\n",
    "                      + \" -r \"+str(res-1), \n",
    "                      shell=True, stdout=subprocess.PIPE)\n",
    "\n",
    "c=subp.stdout.readline()\n",
    "while c:\n",
    "    print(c)\n",
    "    c=subp.stdout.readline()\n",
    "\n",
    "print(subp.returncode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "import open3d as o3d\n",
    "pcd = o3d.io.read_point_cloud(input_file)\n",
    "o3d.visualization.draw_geometries([pcd])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcd_rec = o3d.io.read_point_cloud(output_file)\n",
    "o3d.visualization.draw_geometries([pcd_rec])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## select optimal rho"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "from myutils.pc_error_wrapper import pc_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_rho(item, input_file, output_file, input_file_n, \n",
    "               cubes_d, points_numbers_d, cube_positions_d, scale, cube_size, res):\n",
    "    steps = [0.05]*6+[0.10]*3\n",
    "    MAX = 0\n",
    "    rho = 1.0\n",
    "    optimal_rho = 1.0\n",
    "\n",
    "    for i, step in enumerate(steps):\n",
    "        print('===== select rho =====')\n",
    "        print(rho)\n",
    "        postprocess(output_file, cubes_d.numpy(), \n",
    "                    points_numbers_d, cube_positions_d, scale, cube_size, rho)\n",
    "        results = pc_error(input_file, output_file, input_file_n, res, show=False)\n",
    "        \n",
    "        # record results.\n",
    "        results[\"n_points\"] = get_points_number(output_file)\n",
    "        results[\"rho\"] = rho\n",
    "        if i == 0:\n",
    "            all_results = results.copy(deep=True)\n",
    "        else:\n",
    "            all_results = all_results.append(results, ignore_index=True)\n",
    "        \n",
    "        PSNR = float(results[item])\n",
    "        print('===== results: ', i, rho, item, PSNR)\n",
    "        MAX = max(PSNR, MAX)\n",
    "        if PSNR < MAX:\n",
    "            break\n",
    "        else:\n",
    "            optimal_rho = rho\n",
    "        if item == \"mseF,PSNR (p2point)\":\n",
    "            rho += step\n",
    "        elif item == \"mseF,PSNR (p2plane)\":\n",
    "            rho -= step\n",
    "        else:\n",
    "            print('ERROR', item)\n",
    "            break\n",
    "    return optimal_rho, all_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== select rho =====\n",
      "1.0\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.4908s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.2203s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 21.9627\n",
      "===== results:  0 1.0 mseF,PSNR (p2point) 67.7148\n",
      "===== select rho =====\n",
      "1.05\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3364s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.2875s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 22.3718\n",
      "===== results:  1 1.05 mseF,PSNR (p2point) 68.1193\n",
      "===== select rho =====\n",
      "1.1\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3391s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.4846s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 22.9781\n",
      "===== results:  2 1.1 mseF,PSNR (p2point) 68.5022\n",
      "===== select rho =====\n",
      "1.1500000000000001\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3451s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.5812s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 24.0211\n",
      "===== results:  3 1.1500000000000001 mseF,PSNR (p2point) 68.5065\n",
      "===== select rho =====\n",
      "1.2000000000000002\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3569s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.8998s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 24.8581\n",
      "===== results:  4 1.2000000000000002 mseF,PSNR (p2point) 68.327\n"
     ]
    }
   ],
   "source": [
    "optimal_d1, all_results = select_rho(\"mseF,PSNR (p2point)\", input_file, output_file, input_file_n, \n",
    "                        cubes_d, points_numbers_d, cube_positions_d, scale, cube_size, res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0GElEQVR4nO3dd5wV1f3/8dd7WXpZOiLSRRBR2lpAxIKV2GMsKJbEIGoIosZo/JroL5oYA0bsEkuCAiKo0UTFXomgLCxIb9K7VBFYdvfz+2NmzXW9y1J27mz5PB8PHnvnzJyZz517mc+dc2bmyMxwzjnnCkuLOwDnnHOlkycI55xzSXmCcM45l5QnCOecc0l5gnDOOZeUJwjnnHNJeYJwkZJ0gqR5Ea37H5LujWr5YtbVXlK2pG2Sfl0S64yapLslvRC+biHpW0mVYorlSUl3xbFtt/c8QbgSJckkHVowbWafmln7OGOKyG3Ah2ZW28wePtCVhQfv3WHC2SZpvqRHJTVNWKaKpPGSloT7+aT93Z6ZLTOzWmaWV0xcV0v6bH+3U9Q6zGygmf3xQNbroucJwlUIktJLeJUtgVklHMtYM6sN1AcuAA4CshKTBPAZcAWwZn+2nWoR7HeXQp4g3I9IOlzSR5I2S5ol6dyEef8ImwfeDX/pfiypZTjvk3Cx6WHzxSWSTpK0IqH+Ekm/kTRD0nZJz0hqIumtcH3vSaqXsPw4SWskbZH0iaQj9vI9XC1poqS/SfoGuDucVU/SG+G2Jktqm1Cnp6Qvw219KalnEev+ADgZeDR8n4dJypA0UtJ6SUsl/Z+ktGJiScrMdpvZLOASYD1wS1ieY2YPmdlnwB5/+YfbbR1+PtskvQs0TJjXKjwLSU+IcXG47NeSLpd0OPAk0CN8n5vDZfflvY4tYh0/aO6T9EtJCyVtlPS6pIMT5pmkgZIWhN/JxySpuPfvDpwnCPcDkioD/wbeARoDg4BRkhKbiS4H/khwwMkGRgGYWe9wfuew+WJsEZv5KXAacBhwDvAW8DugEcF3MrFN/y2gXRjL1IJt7aVjgcVAE+C+sOxS4B6gHrCwoFxSfeAN4GGgAfAg8IakBoVXamanAJ8Cvwrf53zgESADaAOcCFwJXFNMLHsUNv+8Bpyw1+/4h0YDWQSf0x+Bq5ItJKkmwfs+KzyD6Qlkm9kcYCDwefg+64ZV9uW9XlHEOhK3fwrwZ+BioCmwFHix0GJnA0cDR4XLnbG3O8HtP08QrrDjgFrA/eEv1g+A/wCXJSzzhpl9Yma7gDsJfh0234dtPGJma81sJcGBdrKZTTOzncCrQNeCBc3sWTPbFm7rbqCzpIy93M4qM3vEzHLNbEdY9qqZfWFmuQTJpktY/hNggZk9Hy4/BphLkMD2SEFH76XAHWGsS4BhQP9iYtmr90DQ5LRPJLUgOKDeZWa7zOwTgsRflHygk6TqZrY6PINJtt4o3uvlwLNmNjX8nO8g+E61SljmfjPbbGbLgA/53+fmIuQJwhV2MLDczPITypYCzRKmlxe8MLNvgY1hvb21NuH1jiTTtSA4GEm6X9IiSVuBJeEyDdk7y5OUJbbdf1ewLYL4lxZatvD7LkpDoHKh+kXus33UjGD/7quDgU1mtr1QTD8SLnMJwS/91WETXIci1hvFe/3Bvg+/U98UWmdRn5uLkCcIV9gqoHlBm3KoBbAyYfr7swVJtQh+4a6KIJZ+wHnAqQRNGq0KNruX9fflUcWrCDqeExV+30XZAOwuVL9w3X1+bHL4GZxDcJa1r1YT9LfULBRTUmb2tpmdRtDEMxf4e8GsQovuz3st7r3/YN+HMTdg7/a9i5AnCFfYZIJfaLdJqqzgUspz+GGbcF9JvSRVIWjbnmRmBb8a1xK0TZeE2sAugl+TNYA/ldB6k3kTOExSP0npki4BOhI0r+1R2FfwEnCfpNoKOu1vBl7Yn0DC7R8OjCG4kunBhHlVJVULJ6tIqpasw9bMlgJTgHsUXB7biyKay8KLBM4LD8y7gG8Jmpwg+DwPCT/r/X2vP1hHEmOAayR1kVSV4HOeHDZfuRh5gnA/YGY5BAeSswh+LT4OXGlmcxMWGw38gaDpoztBR2SBu4F/hlebXHyA4YwkaHpYCcwGJh3g+opkZt8QdITeQpCQbgPONrMNe7mKQcB2gs7Zzwj20bP7GMYlkr4FtgCvh3F0N7PEs7N5BM1wzYC3w9eFz3wK9CPoMN5I8HmNLGK5NIKD/Kpw2ROB68N5HxBczrtGUsG+2Nf3mmwd3zOz94C7gJcJznzaEvRzuJjJBwxy+0LSP4AVZvZ/ccfinIuWn0E455xLyhOEc865pLyJyTnnXFJ+BuGccy6pcvUgrYYNG1qrVq3iDsM558qMrKysDWbWKNm8cpUgWrVqxZQpU+IOwznnygxJSe+wB29ics45VwRPEM4555LyBOGccy4pTxDOOeeS8gThnHMuKU8QzjnnkvIE4ZxzLqlI74OQVBd4GuhEMGjIzwkeT/wkUA3IBW4wsy+S1M0Dvgonl5nZuVHG6lxFlp9v5OYb+Rb8zcsz8szIzc8nL99+9C83cdrCsrz/1c//fpl88vL58XqS1Elcd7vGtfjJkU1JS9vbsaFcFKK+UW44MMHMLgoHC6lBMNjIPWb2lqS+wAPASUnq7jCzLhHH51zkZq7cwqTF33x/8MtPPMD+6EAZHFDz8vN/sGy+BcsUrlPwOvHgGtTJJ9/CA3PCcj84sCfUKY2enfg195x7BEcdUjfuUCqsyBJEOLB8b+Bq+H4gmhxJBtQJF8sgmqEqnSsV5qzeys+e/Jwdu/N+NC9NkJ6WRlpa+FeQXimNSmmikhT8TRPpaf97XTCdVvBXolrlNCqlpVFJBH/D9e2pTnqaqFQp2E7ivIL6ldLSflgexpNe6X/105LG9r/6yeokiyexfiUFy74ybSX3vzWX8x6byCWZzfnNGe1pUKtqDJ9gxRbZ01wldQFGEIwE1hnIAgYTjF/7NsG4wmlAz3B4xML1c4Fsgmao+83sX0VsZwAwAKBFixbdly4t8q5x51Jq83c5nPvoRHbl5jHuup40rF3lBwfKJCOFugTbdu7m4fcX8NzEJVSvUombTzuM/se1JL2Sd52WJElZZpaZdF6ECSKTYIjI481ssqThwFaCs4aPzezlcEjKAWZ2apL6zcxspaQ2BEMW9jGzRXvaZmZmpvmzmFxpkJdv/PwfX/LfRRt4cUAPuresF3dIZdbCddu459+z+XTBBto3qc3d5x5Bj7YN4g6r3NhTgogyFa8gGJpycjg9HugGXAW8EpaNA45JVtnMVoZ/FwMfAV0jjNW5EvXgu/P4eP567jm3kyeHA3Ro49qM/PkxPNW/O9tzcrns75O4cfRUVm3eEXdo5V5kCcLM1gDLJbUPi/oQNDetIhgUHeAUYEHhupLqSaoavm4IHB/Wda7UmzBzNY99uIjLjmlOv2NbxB1OuSCJM444iPduPpEhpx7Ge7PXcsqwj3jk/QXsTNK/40pGpCPKhf0QTwNVgMXANcARBFc3pQM7CS5zzQqbpAaa2bWSegJPAfkESewhM3umuO15E5OL24K12zj/sYm0a1KbsdcdR9X0SnGHVC6t2PQd970xh7dmrqFF/RrcdXZHTj28sffr7IdY+iDi4AnCxWnrzt2c9+hEtu3M5T+DenFQRrW4Qyr3Ji7cwB9en8XCdd9y4mGN+MM5HWnTqFbcYZUpcfVBOFdh5OcbQ17MZvnG73j88m6eHFLk+EMb8tbgE7jr7I5MXbqJMx76hD+/NYdvd+XGHVq54AnCuRIw/P0FvD93Hb8/pyPHtK4fdzgVSuVKafyiV2s+uPUkzu/SjKc+XswpQz/i1WkrKE8tJHHwBOHcAXp39lqGv7+Ai7ofQv/jWsYdToXVqHZV/vqzzrx6Q0+aZlRjyNjp/OzJz5m5ckvcoZVZniCcOwCL1n/LzWOzObJZBvee38k7SUuBri3q8eoNx/PAT4/i6w3bOefRz7jz1a/YtD0n7tDKHE8Qzu2nbTt3M2DkFCqnp/Fk/+5Uq+xXLJUWaWni4qOb88GtJ3F1z1a8+OVyThr6Ec9/voTcvPy4wyszPEE4tx/y841bx01nyTff8Wi/rjSrWz3ukFwSGdUr84dzjuDNX59Ax6Z1uOu1WZzz6ES++Hpj3KGVCZ4gnNsPT3y8iLdnreV3fQ+nZ9uGcYfjitH+oNqM/uWxPH55N7bu2M3FT33Or8dMY82WnXGHVqp5gnBuH304bx1D35nHeV0O5ufHt4o7HLeXJNH3yKa8d/OJ/PqUQ5kwaw2nDPuIxz9ayK5cvxs7GU8Qzu2DJRu2M3jMNDocVIf7LzzKO6XLoOpVKnHz6e15b8iJ9Dq0IQ9MmMeZD33Kh3PXxR1aqeMJwrm9tH1XLtc9n0VamhjRvzvVq3indFnWokENRlyZycifH4ME1/zjS37+jy9ZsmF73KGVGp4gnNsLZsZtL89gwbptPHJZV5rXrxF3SK6E9D6sERMG9+Z3fTswefE3nP63T3hgwly2+93YniCc2xsjPlnMGzNWc9uZHTihXaO4w3ElrEp6GgN6t+XDW0/i7M5NefyjRfQZ9jGvZa+s0Hdje4JwrhifLljPXybM5SdHNuW63m3iDsdFqHGdajx4cRdevr4HDWtXYfCL2VwyYhKzV22NO7RYeIJwbg+Wb/yOQWOm0a5xbR64yDulK4ruLevz2o29+NMFR7Jg7TbOfuRTfv/aTDZ/V7HuxvYE4VwRduTkcd3zWeTnG0/1707Nqulxh+RSqFKa6HdsCz669WT6H9eSFyYt5eShHzFq8lLy8itGs5MnCOeSMDPueGUGc9ZsZfhlXWnVsGbcIbmYZNSozD3ndeKNX59Auya1ufPVmZz32GdkLS3/d2NHmiAk1ZU0XtJcSXMk9ZDURdIkSdmSpkhKOiZ1WL+OpBWSHo0yTucKe3biEv6VvYpbTjuMk9s3jjscVwoc3rQOYwccx8OXdWXDthx++sTnDBmbzbqt5fdu7KjPIIYDE8ysA9AZmAM8ANxjZl2A34fTRfkj8EnEMTr3A/9dtIE/vTmHM45owg0nHRp3OK4UkcS5nQ/mg1tP5MaT2/LGjNWcPPQjnvp4ETm55e8hgJElCEkZQG/gGQAzyzGzzYABdcLFMoBVRdTvDjQB3okqRucKW7l5B4NGT6N1w5oMu7gLaWneKe1+rEaVdH5zRgfeGdKb49o04M9vzeXM4Z/w8fz1cYdWoqI8g2gNrAeekzRN0tOSagI3AX+VtBwYCtxRuKKkNGAYcGtxG5E0IGyqmrJ+ffn6cFxq7dydx8Dns8jJzeep/t2p5Z3SrhitGtbkmauP5rmrj8YMrnr2C345cgrLvvku7tBKRJQJIh3oBjxhZl2B7cDtwPXAEDNrDgwhPMMo5AbgTTNbUdxGzGyEmWWaWWajRn4Dk9s/Zsadr87kq5VbePCSLrT1ge/dPji5Q2Mm3HQCvz2zAxMXbuDUv33MsHfmsSOnbD8EUFHdJSjpIGCSmbUKp08gSBC9gLpmZgouKt9iZnUK1R0FnADkA7WAKsDjZnb7nraZmZlpU6ZMKfH34sq/kZ8v4fevzWJwn3YMOe2wuMNxZdiaLTv581tzeC17FQdnVOPOn3Sk75EHldp7aCRlmVlmsnmRnUGY2RpguaT2YVEfYDZBn8OJYdkpwIIkdS83sxZhcrkVGFlccnBuf325ZCP/79+z6dOhMYP7tIs7HFfGHZRRjeGXduWl63qQUaMKN46eSr+/T2bemm1xh7bPor6KaRAwStIMoAvwJ+CXwDBJ08PpAQCSMiU9HXE8zv3Ami07uf6FqTSvX4O/Xeqd0q7kHNO6Pv8Z1Is/nt+JOWu20vfhT7n79Vls2bE77tD2WmRNTHHwJia3L3bl5nHJU5NYsHYb/7rxeNo1qR13SK6c2rQ9h2HvzmP05GXUrVGF285oz8WZzUvFD5JYmpicK+3ufn0W2cs3M+zizp4cXKTq1azCvecfyeu/6kWbhjW5/ZWvOP/xiUxbtinu0PbIE4SrkEZPXsaYL5Zz48ltObNT07jDcRVEp2YZjBvYg4cu6cKaLTu54PH/cuu46azftivu0JLyBOEqnKylm/jD6zPpfVgjbj6tffEVnCtBkji/azM+uPUkrjuxDa9lr+SUoR/x9KeL2Z1Xuu7G9gThKpR1W3dy/QtZNM2ozsOXdqFSKWgDdhVTrarp3HHW4Uy4qTfdWtbj3jfmcNbwT/lswYa4Q/ueJwhXYeTk5nPDqKls25nLU/27U7dGlbhDco62jWrxj2uO5ukrM8nJzeeKZyYz8Pkslm+M/25sf5aAqzDufWM2U5Zu4pHLunJ40zrFV3AuRSRxascm9GrXkKc/XcxjHy7iw3nruP6ktgw8sS3VKleKJS4/g3AVwktTljPy86UM6N2GczofHHc4ziVVrXIlfnVKO96/5URO69iEh95bQJ9hHzNh5upYxsb2BOHKvenLN/N//5rJ8Yc24LYzvFPalX4H163Oo/26MeaXx1GrajoDX5hK/2e+YOG61N6N7QnClWsbvt3FwBeyaFSrKo9c1o30Sv6Vd2VHj7YNeOPXvbj7nI7MWLGZMx/6lHv/M5utO1NzN7b/b3Hl1u68fG4cNZWN23N4qn936tf0TmlX9qRXSuPq41vz4a0n8bPMQ3hm4tecMvRjxk1ZTn7EY2N7gnDl1p/fnMvkrzdy/0+PpFOzjLjDce6ANKhVlT9feBSv3Xg8zetX5zfjZ3DhE/9l+vLNkW3TE4Qrl16dtoJnJ37NNce34oKuh8QdjnMl5qhD6vLywJ4M/VlnVmzawfmPT+S342ewc3fJjz3hl7m6cmfmyi3c/vJXHNu6Pr/re3jc4ThX4tLSxEXdD+GMI5rw8PsL+GrlFqqml/zvfU8QrlzZuD2H657Pon7NKjx2eTcqe6e0K8dqV6vMnT/pSF6+RTIgkScIV27k5uUzaMxU1n+7i3HX9aBhrapxh+RcSkT1yBj/eeXKjb++PY+JC7/h3vM70bl53bjDca7MizRBSKorabykuZLmSOohqYukSZKyJU2RdEySei0lTQ2XmSVpYJRxurLv39NX8dQni+l/XEsuzmwedzjOlQtRNzENByaY2UWSqgA1gJeAe8zsLUl9gQeAkwrVWw30MLNdkmoBMyW9bmarIo7XlUFz12zltvEzyGxZj7vO7hh3OM6VG5ElCEkZQG/gagAzywFyJBlQ8KS0DOBHB/1w2QJV8aYwV4TN3+UwYGQWtaul8/gV3agSwZUczlVUUZ5BtAbWA89J6gxkAYOBm4C3JQ0lOPD3TFZZUnPgDeBQ4DdFnT1IGgAMAGjRokUJvwVXmuXlG4NfzGb1lh28OKAHjWtXizsk58qVKH9upQPdgCfMrCuwHbgduB4YYmbNgSHAM8kqm9lyMzuKIEFcJalJEcuNMLNMM8ts1KhRFO/DlVIPvjuPj+ev555zO9G9Zb24w3Gu3IkyQawAVpjZ5HB6PEHCuAp4JSwbB/yokzpReOYwEzghojhdGTRh5moe+3ARlx3TnH7H+pmjc1GILEGY2RpguaSC5yv3AWYT9DmcGJadAiwoXFfSIZKqh6/rAb2AeVHF6sqWBWu3cctL0+nSvC53n3tE3OE4V25FfRXTIGBUeAXTYuAa4DVguKR0YCdh/4GkTGCgmV0LHA4MCzu0BQw1s68ijtWVAVt37mbA81lUr5LOk1d0p2p6PCNtOVcRRJogzCwbyCxU/BnQPcmyU4Brw9fvAkdFGZsre/LzjSEvZrN843eM/uVxHJThndLORcmvCXRlxsMfLOD9uev4/TkdOaZ1/bjDca7c8wThyoR3Z6/lofcW8NNuh9D/uJZxh+NcheAJwpV6i9Z/y81jszmyWQb3XdApkqdWOud+zBOEK9W+3ZXLdc9nUTk9jSf7d6daZe+Udi5V/HHfrtTKzzdueSmbrzds5/lfHEOzutXjDsm5CsXPIFyp9cTHi3h71lruOKsDPds2jDsc5yocTxCuVPpw3jqGvjOP87oczC96tY47HOcqJE8QrtRZ+s12Bo+ZRoeD6nD/hUd5p7RzMfEE4UqV7btyGTAyi7Q0MaJ/d6pX8U5p5+LiCcKVGmbGbS/PYMG6bTxyWVea168Rd0jOVWieIFypMeKTxbwxYzW3ndmBE9r5o9udi5snCFcqfLZgA3+ZMJefHNmU63q3iTsc5xyeIFwpsHzjd/xqzFTaNa7NAxd5p7RzpYUnCBerHTl5XPd8Fvn5xlP9u1Ozqt+76Vxp4f8bXWzMjDtemcGcNVt59qqjadWwZtwhOecS+BmEi81zE5fwr+xV3HLaYZzcoXHc4TjnCok0QUiqK2m8pLmS5kjqIamLpEmSsiVNkfSjManDZT6XNEvSDEmXRBmnS73PF33DfW/O4YwjmnDDSYfGHY5zLomom5iGAxPM7KJw2NEawEvAPWb2lqS+wAPASYXqfQdcaWYLJB0MZEl628w2RxyvS4GVm3fwq9FTadWgBkN/1pm0NO+Udq40iixBSMoAegNXA5hZDpATjjNdJ1wsA1hVuK6ZzU94vUrSOqARsDmqeF1q7Nydx/UvZLErN58RV2ZSu1rluENyzhUhyjOI1sB64DlJnYEsYDBwE/C2pKEETVw997SSsAmqCrCoiPkDgAEALVq0KKnYXQTMjP/710xmrNjC36/MpG2jWnGH5Jzbgyj7INKBbsATZtYV2A7cDlwPDDGz5sAQ4JmiViCpKfA8cI2Z5SdbxsxGmFmmmWU2auR335ZmL0xayvisFQzu047TOjaJOxznXDGiTBArgBVmNjmcHk+QMK4CXgnLxgE/6qQGkFQHeAO408wmRRinS4Evl2zknn/Ppk+Hxgzu0y7ucJxzeyGyBGFma4DlktqHRX2A2QR9DieGZacACwrXDTu0XwVGmtn4qGJ0qbFmy06uf2EqzevX4MFLunintHNlRNRXMQ0CRoUH/MXANcBrwHBJ6cBOwv4DSZnAQDO7FriYoIO7gaSrw3VdbWbZEcfrStiu3DyuH5XFjpxcxvzyWDKqe6e0c2VFpAkiPKBnFir+DOieZNkpwLXh6xeAF6KMzaXG3a/PZtqyzTx5RTfaNakddzjOuX3gd1K7yIyevIwxXyzjxpPbcmanpnGH45zbR54gXCSmLtvEH16fSe/DGnHzae2Lr+CcK3U8QbgSt27bTq5/IYumGdV5+NIuVPJOaefKJH+aqytRObn53DhqKlt35PLKDcdQt0aVuENyzu0nTxCuRN37xmy+XLKJRy7ryuFN6xRfwTlXau13E5Mkf3i/+4FxU5Yz8vOlDOjdhnM6Hxx3OM65A1RsgpDUTFJmeC8DkhpL+hNJbnBzFdeMFZu5818zOf7QBtx2hndKO1ce7DFBSLoJyAYeASZJuhaYA1Qnyb0MrmLa8O0uBj6fRaNaVXnksm6kV/JrH5wrD4rrgxgAtDezjZJaAPOB480sK/rQXFmwOy/olP5mew4vX9+T+jW9U9q58qK4n3o7zWwjgJktA+Z5cnCJ/vzmXCZ/vZH7f3oknZplxB2Oc64EFXcGcYikhxOmmyZOm9mvownLlQX/mraSZyd+zTXHt+KCrofEHY5zroQVlyB+U2jazx4cADNXbuH2V2ZwbOv6/K7v4XGH45yLwB4ThJn9M1WBuLJj4/Ycrns+i3o1qvDY5d2o7J3SzpVLe0wQkp4DrIjZZma/KPmQXGmWl28MfnEa67/dxbjretCwVtW4Q3LORaS4Jqb/JCkrGCq0UsmH40q7Rz9YyKcLNnD/hUfSuXnduMNxzkWouCamlwteS2oD/I5gIJ/72cNY0q58+u/CDTz0/nwu7NqMS45uHnc4zrmI7c2d1B0kvQD8m2Cwn45m9oSZ5exF3bqSxkuaK2mOpB6SukiaJClb0hRJRY1JPUHSZknJzmJciq3btpNfv5hN20a1uPeCTkj+hFbnyrvi+iDGEdwxPYygWSkPqFNwcCi4R2IPhgMTzOyi8FEdNYCXgHvM7C1JfYEHgJOS1P1ruPx1e/1uXCTy8o3BY7LZviuX0b88lhpV/BmPzlUExf1PP5qgk/pW4BYg8WejAW2Kqigpg6A56mqA8IwjR5IBBY/5zABWJatvZu9LOqm4N+CiN/y9+Xy++BuG/qwzh/mwoc5VGMX1QbQ6gHW3BtYDz0nqTHAPxWDgJuBtSUMJmrh6HsA2kDSA4JEgtGjR4kBW5ZL4ZP56HvlwIT/rfggXdfeb4ZyrSIp7WF/L8EygYPpkScMlDSl4uusepAPdgCfMrCuwHbgduB4YYmYFV0MdUGe3mY0ws0wzy2zUqNGBrMoVsnbrToaMzeawxrX5f+d1ijsc51yKFddJ/RJQE0BSF2AcsAzoAjxeTN0VwAozmxxOjydIGFcBr4Rl44CkndQuXrl5+QwaPY0du/N47PJuVK/iVzU7V9EUlyCqm1lBH8EVwLNmNgy4hmIO7Ga2BlguqWBwgD7AbII+hxPDslPwcSVKpb+9N58vlmzkvgs6cWjjWnGH45yLQXGd1Imd0qcAdwCYWf5eXuY4CBgVNkctJkgsrwHDJaUDOwn7DyRlAgPN7Npw+lOgA1BL0grgF2b29t6+Mbf/Ppq3jsc+XMSlRzf3h/A5V4EVlyA+kPQSsBqoB3wAIKkpUOx9EGaWDWQWKv6MJIMNmdkU4NqE6ROKW78reau37GDI2Gw6HFSbu889Iu5wnHMxKi5B3ARcAjQFepnZ7rD8IODOCONyMdgd9jvk5Obz+OXdqFbZ+x2cq8iKu8zVgBeTzJoOXBZJRC42Q9+Zx5Slm3j4sq60aeT9Ds5VdMVd5lpH0h2SHpV0ugKDCPoTLk5NiC4V3p+zlqc+Xszlx7bg3M4Hxx2Oc64UKK6J6XlgE/A5Qf/A7wg6rs8P+xdcObBy8w5uGTedjk3rcNfZHeMOxzlXShSXINqY2ZEAkp4m6KxuYWY7I4/MpURObj6/Gj2V3Dzzfgfn3A8Udx9EQac0ZpZHcOObJ4dy5K9vz2Xass385adH0aphzbjDcc6VIsWdQXSWtDV8LaB6OC2CPuw6RVd1pd27s9fy90+/5soeLfnJUU3jDsc5V8oUdxWTtzeUU8s3fsctL2VzZLMM7vzJ4XGH45wrhXy0+QooJzefX42ZhgGP9etG1XT/HeCc+zEf+aUC+vNbc5i+fDNPXtGdFg1qxB2Oc66U8jOICmbCzNU8N3EJ1xzfijM7HRR3OM65UswTRAWy7Jvv+M34GXQ+JIM7zvJ+B+fcnnmCqCB25eZx4+ipCHi0XzeqpPtH75zbM++DqCD+9MYcvlq5hRH9u9O8vvc7OOeK5z8jK4A3Zqzmn58v5dperTn9CO93cM7tHU8Q5dySDdv57csz6NqiLr89q0Pc4TjnypBIE4SkupLGS5oraY6kHpK6SJokKVvSFElJhy6VdJWkBeG/q6KMs7zauTuPG0ZNJb2SeLRfNypX8t8Dzrm9F3UfxHBggpldFA47WgN4CbjHzN6S1Bd4ADgpsZKk+sAfCEajMyBL0utmtinieMuVP/5nNrNXb+XZqzNpVrd63OE458qYyH5SSsoAegPPAJhZjpltJjjgFzzDKQNYlaT6GcC7ZrYxTArvAmdGFWt59Fr2SkZNXsZ1J7bhlA5N4g7HOVcGRXkG0RpYDzwnqTOQBQwmGMb0bUlDCRJUzyR1mwHLE6ZXhGU/ImkAMACgRYsWJRV7mbZo/bf87pWvyGxZj1tPbx93OM65MirKRul0oBvwhJl1BbYDtwPXA0PMrDkwhPAMY3+Z2QgzyzSzzEaNGh1ozGXezt153DhqKlXS03ikX1fvd3DO7bcojx4rCMaPmBxOjydIGFcBr4Rl44BkndQrgeYJ04eEZa4Yd78+i7lrtvHgJV1omuH9Ds65/RdZgjCzNcBySQVtHH2A2QR9DieGZacAC5JUfxs4XVI9SfWA08MytwevTlvBi18u54aT2nJy+8Zxh+OcK+OivoppEDAqvIJpMXAN8BowXFI6sJOw/0BSJjDQzK41s42S/gh8Ga7n/5nZxohjLdMWrvuWO1+dyTGt63PzaYfFHY5zrhyQmcUdQ4nJzMy0KVOmxB1Gyu3IyeP8xyay4dtdvDn4BJrUqRZ3SM65MkJSlpllJpvnz2IqB37/2kzmr9vGyJ8f48nBOVdi/BKXMm581grGZa1g0MmHckI7v4rLOVdyPEGUYfPXbuP//vUVx7Wpz+BTvd/BOVeyPEGUUdt35XLDqKnUqlqZhy/tSqU0xR2Sc66c8T6IMsjMuOtfM1m0/lte+MWxNPZ+B+dcBPwMogwaN2UFr0xbyeA+7Tj+0IZxh+OcK6c8QZQxc9ds5a7XZtLr0IYMOqVd3OE458oxTxBlyLdhv0Od6pX52yVdvN/BORcpTxBlhJlx56tfsWTDdh65rCuNaleNOyTnXDnnCaKMGPPFcl7LXsXNpx3GcW0axB2Oc64C8ARRBsxatYW7/z2LE9o15IaTDo07HOdcBeEJopTbtnM3N46aSr0alXnoki6keb+Dcy5F/D6IUszMuP2Vr1i+aQdjfnkcDWp5v4NzLnX8DKIUe2HSUt6YsZpbTj+MY1rXjzsc51wF4wmilJq5cgt//M8cTm7fiIG928YdjnOuAvIEUQpt3bmbG0ZNpUGtKgy72PsdnHPxiDRBSKorabykuZLmSOohaayk7PDfEknZRdQdLGmmpFmSbooyztLEzPjt+Bms2ryDR/t1pX7NKnGH5JyroKLupB4OTDCzi8JhR2uY2SUFMyUNA7YUriSpE/BL4BggB5gg6T9mtjDieGP3z/8u4a2Za/hd3w50b+n9Ds65+ER2BiEpA+gNPANgZjlmtjlhvoCLgTFJqh8OTDaz78wsF/gYuDCqWEuL6cs3c9+bc+jToTHX9moTdzjOuQouyiam1sB64DlJ0yQ9LalmwvwTgLVmtiBJ3ZnACZIaSKoB9AWaJ9uIpAGSpkiasn79+pJ+Dymz5bvd3Dh6Ko1rV2PYxZ2938E5F7soE0Q60A14wsy6AtuB2xPmX0byswfMbA7wF+AdYAKQDeQVsewIM8s0s8xGjcrmkJtmxq3jp7Nmy04e6deVujW838E5F78oE8QKYIWZTQ6nxxMkDCSlEzQZjS2qspk9Y2bdzaw3sAmYH2GssXrms695d/Zabj+rA91a1Is7HOecAyJMEGa2BlguqX1Y1AeYHb4+FZhrZiuKqi+pcfi3BUEyGR1VrHGatmwT9781l9M7NuEXvVrHHY5zzn0v6quYBgGjwiuYFgPXhOWXUqh5SdLBwNNm1jcsellSA2A3cGNiB3d5sfm7HH41ehoHZVTjrxd1Jui3d8650iHSBGFm2UBmkvKrk5StIuiMLpg+IcrY4mZm3DpuOuu27WT8wJ5k1Kgcd0jOOfcDfid1TP7+6WLem7OOO/seTufmdeMOxznnfsQTRAyylm7kLxPmcVang7iqZ6u4w3HOuaQ8QaTYxu1Bv0OzutX5y0VHeb+Dc67U8vEgUig/37j5pWy++TaHV27oSZ1q3u/gnCu9/AwihZ78ZBEfzVvPXWcfTqdmGXGH45xze+QJIkW++Hojw96Zz9lHNeWK41rGHY5zzhXLE0QKfPPtLgaNmUqL+jX484VHer+Dc65M8AQRsfx8Y8hL09n03W4e7deV2t7v4JwrIzxBROzxjxbyyfz13H3OERxxsPc7OOfKDk8QEfp80Tc8+O58zu18MJcdk/Rp5c45V2p5gojI+m27+PWL02jVoCZ/8n4H51wZ5PdBRCAv37hp7DS27tjNyJ8fQ62qvpudc2WPH7ki8MgHC5i48Bv+8tMjObxpnbjDcc65/eJNTCVs4sINDH9/ARd2bcbFmd7v4JwruzxBlKB123Yy+MVs2jaqxb0XdPJ+B+dcmRZpgpBUV9J4SXMlzZHUQ9JYSdnhvyWSsouoO0TSLEkzJY2RVC3KWA9UXr4xeEw223fl8vjl3ahRxVvvnHNlW9RHseHABDO7KBxVroaZXVIwU9IwYEvhSpKaAb8GOprZDkkvEYxC94+I491vw9+bz+eLv2HozzpzWJPacYfjnHMHLLIEISkD6A1cDWBmOUBOwnwBFwOn7CG26pJ2AzWAVVHFeqA+mb+eRz5cyEXdD+Gi7ofEHY5zzpWIKJuYWgPrgeckTZP0tKSaCfNPANaa2YLCFc1sJTAUWAasBraY2TsRxrrf1m7dyZCx2bRrXIs/ntcp7nCcc67ERJkg0oFuwBNm1hXYDtyeMP8yYEyyipLqAecRJJmDgZqSrihi2QGSpkiasn79+pKMv1i5efkMGj2NHbvzePzyblSvUiml23fOuShFmSBWACvMbHI4PZ4gYSApHbgQGFtE3VOBr81svZntBl4BeiZb0MxGmFmmmWU2atSoRN9AcR58dz5fLNnIfRd04tDG3u/gnCtfIksQZrYGWC6pfVjUB5gdvj4VmGtmK4qovgw4TlKNsK+iDzAnqlj3x4fz1vH4R4u49OjmXNDV+x2cc+VP1PdBDAJGSZoBdAH+FJZfSqHmJUkHS3oTIDzrGA9MBb4K4xwRcax7bfWWHdw8NpsOB9Xm7nOPiDsc55yLhMws7hhKTGZmpk2ZMiXSbezOy+eyEZOYs3orrw/qRdtGtSLdnnPORUlSlpllJpvnd3Pto6HvzGPK0k0Mv7SLJwfnXLnmj9rYB+/PWctTHy+m37EtOK9Ls7jDcc65SHmC2EsrN+/glnHT6di0Dr8/u2Pc4TjnXOQ8QeyFnNx8fjV6Krl5xuOXd6NaZb/fwTlX/nkfxF54YMJcpi3bzGP9utGqYc3iKzjnXDngZxDFeGfWGp7+7Guu7NGSnxzVNO5wnHMuZTxB7MHyjd9x67jpHNksgzt/cnjc4TjnXEp5gihCQb+DGTzWrxtV073fwTlXsXgfRBH+/NYcpq/YwhOXd6NFgxpxh+OccynnZxBJTJi5mucmLuHqnq0460jvd3DOVUyeIApZ9s13/Gb8DDofksHv+nq/g3Ou4vIEkWBXbh43jp6KgEf7daNKuu8e51zF5X0QCe57Yw5frdzCiP7daV7f+x2ccxWb/0QO/WfGKkZ+vpRre7Xm9CMOijsc55yLnScI4OsN27n95a/o2qIuvz2rQ9zhOOdcqVDhE8TO3XncOGoqldLEo/26UblShd8lzjkHRJwgJNWVNF7SXElzJPWQNFZSdvhviaTsJPXaJyyTLWmrpJuiiDHfjA5Na/PgxZ1pVrd6FJtwzrkyKepO6uHABDO7SFIVoIaZXVIwU9IwYEvhSmY2j2CIUiRVAlYCr0YRYI0q6Tx4cZcoVu2cc2VaZAlCUgbQG7gawMxygJyE+QIuBk4pZlV9gEVmtjSaSJ1zziUTZRNTa2A98JykaZKelpT4rOwTgLVmtqCY9VwKjIkqSOecc8lFmSDSgW7AE2bWFdgO3J4w/zKKOfCHzVLnAuP2sMwASVMkTVm/fv2BR+2ccw6INkGsAFaY2eRwejxBwkBSOnAhMLaYdZwFTDWztUUtYGYjzCzTzDIbNWpUAmE755yDCBOEma0BlktqHxb1AWaHr08F5prZimJWU+xZhnPOuWhEfdH/IGCUpBkEVyX9KSz/Ub+CpIMlvZkwXRM4DXgl4hidc84lEellrmaWDWQmKb86SdkqoG/C9HagQYThOeec2wO/bdg551xSMrO4YygxktYD+3u/RENgQwmGU1I8rn3jce0bj2vflMe4WppZ0it8ylWCOBCSppjZj5rD4uZx7RuPa994XPumosXlTUzOOeeS8gThnHMuKU8Q/zMi7gCK4HHtG49r33hc+6ZCxeV9EM4555LyMwjnnHNJeYJwzjmXVLlMEJKelbRO0swi5kvSw5IWSpohqVvCvKskLQj/XZVQ3l3SV2Gdh8PxLFISl6Qukj6XNCssTxx06R+Svk4Yfa9LquIK5+UlbPv1hPLWkiaHdcaGT+ZNSVySTi40IuFOSeeH81KxvzqEn9cuSbcWmnempHlhzLcnlKdifyWNS1JzSR9Kmh1+xwYnzLtb0sqE/dU32bqjiCuctyT8f5ctaUpCeX1J74b/T9+VVC9VcWkPI16maH9dHn7fv5L0X0mdE+aV7PfLzMrdP4KBiroBM4uY3xd4CxBwHDA5LK8PLA7/1gtf1wvnfREuq7DuWSmM6zCgXfj6YGA1UDec/gdwURz7K5z3bRF1XgIuDV8/CVyfyrgSlqkPbCQYzTBV+6sxcDRwH3BrQnklYBHQBqgCTAc6pnB/FRVXU6Bb+Lo2MD8hrrsTl03l/grnLQEaJqnzAHB7+Pp24C+pjKvQZ7qG4GazVO2vnvzvuHQW/ztOlPj3q1yeQZjZJwQHhaKcB4y0wCSgrqSmwBnAu2a20cw2Ae8CZ4bz6pjZJAv28Ejg/FTFZWbzLRxYyYJnVq0DSuzZ5gewv5KSJIKRAseHRf8khfur0DIXAW+Z2Xf7uv39jcvM1pnZl8DuQrOOARaa2WILRlh8ETgvVfurqLjMbLWZTQ1fbwPmAM32dfslHVcxziPYT5Di/VVIiY94uRdx/Tc8PgFMAg4JX5f496tcJoi90AxYnjC9IizbU/mKJOWpiut7ko4h+HWwKKH4vvCU82+SqqY4rmoKBmyaVNCMQ/CQxc1mlptk+VTFVSDZiIRR76+iFBVvqvZXsSS1AroCkxOKfxXur2f3pynnABnwjqQsSQMSypuY2erw9RqgSYrjKpDs+5XK/fULgrNoiOD7VVETRJkU/jp+HrjGzPLD4juADgSnwvWB36Y4rJYW3OLfD3hIUtsUb79I4f46Eng7oTju/VVqSaoFvAzcZGZbw+IngLYEj+tfDQxLcVi9zKwbQVPKjZJ6F14gPKtP+fX6Sj7iZcr2l6STCRJEZN/hipogVgLNE6YPCcv2VH5IkvJUxYWkOsAbwJ1hcwrwffOAmdku4DmC08yUxWVmBX8XAx8R/Pr8hqC5J73w8qmKK3Qx8KqZfd9EkKL9VZSi4k3V/iqSpMoEyWGUmX0/BouZrTWzvPAHyd9J7f5K/H6tA15N2P7agubE8O+6VMYV+tGIl6naX5KOAp4GzjOzb8LiEv9+VdQE8TpwpQLHAVvC09W3gdMl1QtPDU8H3g7nbZV0XNiedyXwWqriCn+pvErQ3j4+sULCfxIRtCsmvfIhorjqFTTRSGoIHA/MDn/RfUjQ/g9wFSncXwnzfzQiYYr2V1G+BNqFV5RUIWieeD2F+yupcF88A8wxswcLzUvs07mAFO4vSTUl1S54TfD/sWD7rxPsJ0jx/kpQ5PcrFMn+ktSCYCC1/mY2P2FWyX+/9qYnu6z9I/jQVhN0Lq0gOA0bCAwM5wt4jKAd/ysgM6Huz4GF4b9rEsozCT7sRcCjhHehpyIu4IqwTnbCvy7hvA/CZWcCLwC1UhhXz3B6evj3FwnrbENw5ddCglPwqin+HFsR/EpKK7TOVOyvg8LyrcDm8HWdcF5fgquEFhGcDaZyfyWNC+hF0EQzI+H71Tes83y4v2YQHJSbpjCuNuF3azowq9D+agC8DywA3gPqp/hzrEnwyzyj0DpTsb+eBjYlfFZTEuqW6PfLH7XhnHMuqYraxOScc64YniCcc84l5QnCOedcUp4gnHPOJeUJwjnnXFKeIJyLgKRWRT2N07mywhOEc9EQ/v/LlXH+BXauhIRnDfMkjSS4Ea+6pL8rGGPhHUnVw+W6hA83nCHp1RgegOfcXvEE4VzJagc8DhxB8Fycx8zsCII7cX8aLjMS+K2ZHUVw1+0fYojTuWJ5gnCuZC21/z1M8Wszyw5fZwGtJGUQDPb0cVj+T4IBYpwrdTxBOFeytie83pXwOg9Ix7kyxBOEcylkZluATZJOCIv6Ax/voYpzsfFfNM6l3lXAk5JqEIx7fk3M8TiXlD/N1TnnXFLexOSccy4pTxDOOeeS8gThnHMuKU8QzjnnkvIE4ZxzLilPEM4555LyBOGccy6p/w+iJFEXzUoCDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.array(all_results[\"rho\"][:]), np.array(all_results[\"mseF,PSNR (p2point)\"][:]))\n",
    "plt.title(\"optimal rho for D1 distortion\")\n",
    "_ = plt.xlabel('rho')\n",
    "_ = plt.ylabel('RSNR')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===== select rho =====\n",
      "1.0\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.7957s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.2509s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 22.0556\n",
      "===== results:  0 1.0 mseF,PSNR (p2plane) 71.1413\n",
      "===== select rho =====\n",
      "0.95\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3473s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 3.0599s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 21.2694\n",
      "===== results:  1 0.95 mseF,PSNR (p2plane) 71.3808\n",
      "===== select rho =====\n",
      "0.8999999999999999\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.3454s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 2.8367s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 20.7336\n",
      "===== results:  2 0.8999999999999999 mseF,PSNR (p2plane) 71.5715\n",
      "===== select rho =====\n",
      "0.8499999999999999\n",
      "===== Post process =====\n",
      "Classify and extract points: 0.342s\n",
      "Write point cloud to longdress_vox10_1300_rec.ply: 2.7797s\n",
      "===== measure PCC quality using `pc_error` version 0.13.4 20.1709\n",
      "===== results:  3 0.8499999999999999 mseF,PSNR (p2plane) 71.3023\n"
     ]
    }
   ],
   "source": [
    "optimal_d2, all_results = select_rho(\"mseF,PSNR (p2plane)\", input_file, output_file, input_file_n, \n",
    "                        cubes_d, points_numbers_d, cube_positions_d, scale, cube_size, res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu40lEQVR4nO3deXhU9fn+8feTlSRsYRHZERSUfQkornWtu+KCC7KJItpWbbVWu2lrbf35rba1iIoCsihWLHZx16p1BxI2kU1EQED2nQAJ5Pn9MSdtTCeEkJmcSXK/rmsuTmbOOXPPDDPPnGU+j7k7IiIipSWFHUBERBKTCoSIiESlAiEiIlGpQIiISFQqECIiEpUKhIiIRKUCIXFlZqeY2ZI4rfsZM/tNvOYvZ12dzGyume00s1tjsc54M7P7zGxKMN3GzHaZWXJIWZ4ws1+Ecd9y6FQgJKbMzM3s6OK/3f0Dd+8UZqY4uQt4193rufujlV1Z8OFdGBScnWa21MxGm1nzEvOcYGZvmdkWM9toZtNK3l4R7r7K3eu6+4Fycg0zsw8P5z4Otg53H+Xu91dmvRJ/KhBSK5hZSoxX2Rb4PMZZ/uLu9YBGwADgSCCvRBHIBsYC7YL73wlMOJwMVSUOz7tUIRUI+R9mdpyZvWdm28zsczO7uMRtzwS7B94Kvun+28zaBre9H8w2L9h9cZWZfcfMVpdYfoWZ/djM5pvZbjMbZ2bNzOy1YH1vm1l2ifmnmdk6M9tuZu+bWZdDfAzDzOwjM/uDmW0G7gtuyjazV4L7mmFmHUosc6KZzQrua5aZnVjGut8BTgdGB4+zo5k1MLNJwTf7lWb2czNLKidLVO5e6O6fA1cBG4E7gutfc/dp7r7D3fOB0cBJB3kOjgpen51m9hbQpMRt7YKtvZQSGZcH835lZoPM7DjgCaB/8Di3BfNW5LH+pYx1fGt3n5ndaGbLgq2jf5hZixK3uZmNMrMvgv+Tj5mZHew5lNhQgZBvMbNU4J/Am8ARwA+AZ82s5G6iQcD9RD5w5gLPArj7qcHtPYLdF38p424uB84GOgIXAa8BPwWaEvk/WXKf/mvAMUGW2cX3dYiOB5YDzYAHguuuBn5F5Nv4suLrzawR8ArwKNAYeAR4xcwal16pu58BfAB8P3icS4E/Aw2A9sBpwBBgeDlZDirY/fN34JQyZjmVg2/FPAfkEXmd7geGRpvJzLKIPO7zgi2YE4G57r4IGAV8EjzOhsEiFXms15WxjpL3fwbwO2Ag0BxYCTxfarYLgb5A92C+7x7kcUuMqEBIaScAdYEH3b3A3d8BXgauKTHPK+7+vrvvA35G5Nth6wrcx5/dfb27ryHyQTvD3ee4+17gJaBX8YzuPt7ddwb3dR/Qw8waHOL9rHX3P7v7fnffE1z3krvPdPf9RIpNz+D6C4Av3H1yMP9UYDGRAnZQFjnQezVwT5B1BfAwMLicLIf0GIjscip9n92BXwI/LiNTGyIfqL9w933u/j6Rwl+WIqCrmWW4+zfBFky09cbjsQ4Cxrv77OB1vofI/6l2JeZ50N23ufsq4F3++7pJHKlASGktgK/dvajEdSuBliX+/rp4wt13AVuC5Q7V+hLTe6L8XRciH0Zm9qCZfWlmO4AVwTxNODRfR7luXYnp/OL7IpJ/Zal5Sz/usjQBUkstX+ZzVkEtiTy//2GRkwBeA25z9w/KWK4FsNXdd5fK9D+Cea4i8k3/m2AX3LFlrDcej/Vbz33wf2pzqXWW9bpJHKlASGlrgdbF+5QDbYA1Jf7+z9aCmdUl8g13bRyyXAtcApxFZJdGu+K7PcTlKzJU8VoiB35LKv24y7IJKCy1fOllKzxscvAaXERkK6v4urbA28D97j75IIt/Q+R4S1apTFG5+xvufjaRXTyLgafKyH04j7W8x/6t5z7I3JhDe+4ljlQgpLQZRL6h3WVmqWb2HSIfUiX3CZ9vZiebWRqRfdufunvxt8b1RPZNx0I9YB+Rb5OZwG9jtN5oXgU6mtm1ZpZiZlcBnYnsXjuo4FjBC8ADZlYv+BD/ETDlcIIE938cMJXImUyPBNe3BN4BRrv7E+VkWgnkAr8yszQzO5kydpcFJwlcEnww7wN2EdnlBJHXs1XwWh/uY/3WOqKYCgw3s55mlk7kdZ4R7L6SEKlAyLe4ewGRD5LziHxbHAMMcffFJWZ7DriXyK6PPkQORBa7D5gYnG0ysJJxJhHZ9bAGWAh8Wsn1lcndNxM5EHoHkYJ0F3Chu286xFX8ANhN5ODsh0Seo/EVjHGVme0CtgP/CHL0cffirbMbiBTf+4IzgnYF85flWiIHjLcQeb0mlTFfEpEP+bXBvKcBNwe3vUPkQPg6Myt+Lir6WKOt4z/c/W3gF8BfiWz5dCBynENCZmoYJBVhZs8Aq93952FnEZH40haEiIhEpQIhIiJRaReTiIhEpS0IERGJqkYNpNWkSRNv165d2DFERKqNvLy8Te7eNNptNapAtGvXjtzc3LBjiIhUG2YW9Rf2oF1MIiJSBhUIERGJSgVCRESiUoEQEZGoVCBERCQqFQgREYlKBUJERKJSgZCE9PWWfF7MW01RkYaCEQlLjfqhnNQMW3cXMOjpGazaks+7izfw8MAe1ElNDjuWSK2jLQhJKIUHirjl2dms27GXof3b8spn33DtU5+yede+sKOJ1DoqEJJQfvPyQj5ZvpnfDejGry7pyuODevP52h1c9vjHfLnxYM3TRCTWVCAkYTw3YxUTP1nJyFPbc3mfVgCc1605U0eewK69+7lszMfMWL455JQitYcKhCSEGcs388u/L+C0jk35ybnHfuu23m2yeemWk2hcN43B42bytzlrQkopUruoQEjoVm/N5+ZnZ9OmcSaPXtOL5CT7n3naNM5k+s0n0qtNQ27/y1z+/K8vULMrkfhSgZBQ7d63nxsm5lJ4oIinh+TQICO1zHkbZqYxecTxXNarJQ+/tZQfvzifgv1FVZhWpHbRaa4SmqIi585p81i6ficThvejfdO65S6TlpLEwwN70KZxJn98+wvWbtvD49f1OWhhEZHDoy0ICc2j73zBawvW8dPzj+O0jlEbWkVlZtx+VkcevrIHs1Zs4YrHP+brLflxTCpSO6lASChe++wb/vj2F1zRpxUjTj7qsNZxeZ9WTLr+eNbv2MuAMR8z7+ttsQ0pUsupQEiVW7h2Bz96YR692zTkgQFdMfvfg9KHqn+Hxky/5UQy0pK4auwnvPH5uhgmFandVCCkSm3atY8bJ+XSMDOVJwb3IT2l8kNoHH1EPV665SSOPbI+o6bkMe7Dr3SGk0gMqEBIlSnYX8QtU2azadc+xg7O4Yh6dWK27iZ105l64wmc2+VI7n95Iff943P2H9AZTiKVoQIhVcLdufcfC5i5YgsPXdGdbq0axPw+MtKSeeza3ow8tT0TP1nJTZPz2L1vf8zvR6S2UIGQKjH505VMnfk1t3ynA5f0bBm3+0lKMn56/nHcf2lX3l2ygYFPfsL6HXvjdn8iNZkKhMTdx8s28at/LuSs447gznM6Vcl9Dj6hLeOG9WXFpt1c+thHLF63o0ruV6QmUYGQuFq1OZ9bnptN+yZZ/OGqniRFGUYjXk7vdAQvjOpPkTtXPP4J7y/dWGX3LVITqEBI3Ozat58bJs3CHZ4emkO9OlX/a+cuLRrwt++dROtGmQx/ZhZTZ66q8gwi1ZUKhMRFUZFz+/Nz+XLjbsYM6k3bxlmhZWneIINpo/pz8tFNuGf6Zzz42mK1MhU5BCoQEhePvLWUtxet55cXduako5uEHYe66SmMG5rDoOPb8MS/v+QHz89hb+GBsGOJJDQN1icx9895axn97jKu6deaIf3bhh3nP1KSk/jNpV1p2ziT3766mG+27eGpITk0rpsedjSRhKQtCImpz1Zv58cvzqNvu2x+dXHlhtGIBzNj5KkdGFOilelytTIViUoFQmJmw869jJycS+OsdB6/rg9pKYn73+v8kq1MH/+YmV9tCTuSSMJJ3HewVCv79h9g1OQ8tuUXMnZIH5pUg902xa1MG2Wlcd3TM/j7XLUyFSlJBUIqzd352UsLmL1qG48M7EGXFrEfRiNeSrYyve15tTIVKUkFQipt/EcreDFvNbedeQzndWsedpwKUytTkeh0FpNUyvtLN/LAKws5t8uR3HbmMWHHOWzFrUxbN8rkT/9SK1MR0BaEVMLyjbv4/nOz6disHg8P7FGlw2jEg5nxw7PVylSkmAqEHJYdewu5YVIuKclJPDUkh6z0mrMxqlamIhEqEFJhB4qcW6fOYdXmfB4f1JvWjTLDjhRzamUqEscCYWadzGxuicsOM7vdzK40s8/NrMjMcg6y/Aoz+yxYNjdeOaXiHnp9Me8t2civL+nK8e0bhx0nbtTKVGq7uBUId1/i7j3dvSfQB8gHXgIWAJcB7x/Cak4P1lFmIZGqNX32ap58fzlD+rfl2uPbhB0n7opbmX63839bmR7QQH9SS1TVLqYzgS/dfaW7L3L3JVV0vxJDc1Zt5e7pn9G/fWN+cWHnsONUmYy0ZMYM+m8r05GTctXKVGqFqioQVwNTK7iMA2+aWZ6ZjSxrJjMbaWa5Zpa7caMawsTLuu17uWlyHs3qpzNmUG9Sk2vX4Su1MpXaKO7vcjNLAy4GplVw0ZPdvTdwHvA9Mzs12kzuPtbdc9w9p2nTppVMK9HsLTzATZMj35qfHtKX7Ky0sCOFZvAJbRk3NNLKdIBamUoNVxVfA88DZrv7+oos5O5rgn83EDl20S8O2aQc7s7df53P/DXb+ePVveh0ZL2wI4Xu9GMjrUwPqJWp1HBVUSCuoYK7l8wsy8zqFU8D5xA5uC1V7Mn3l/O3uWu585xOnN25WdhxEkZxK9NW2RlqZSo1VlwLRPDhfjYwvcR1A8xsNdAfeMXM3giub2FmrwazNQM+NLN5wEzgFXd/PZ5Z5X+9s3g9/+/1xVzYvTm3fKdD2HESTvMGGbx484n/aWX6/15XK1OpWawmndedk5Pjubn6yUQsLNuwk0sf+5h2TTKZdtOJZKQlhx0pYe0/UMS9//icZ2es4oLuzXn4yh7USdXzJdWDmeWV9VOCmjM+gsTMtvwCbpiYS53UZMYOzlFxKEfpVqbrtu9l7OA+amUq1V7tOldRyrX/QBHff24Oa7ft5cnBfWjRMCPsSNVCyVamC9ZsVytTqRFUIORbHnh1ER8u28RvBnSlT9vssONUO2plKjWJCoT8x19mrWLCRysYcfJRDMxpHXacakutTKWmUIEQAHJXbOHnf1vAKcc04Z7zjg07TrWnVqZSE6hACGu27WHUlDxaZWcy+prepNSyYTTipWFmGpNG9PtPK9O71MpUqhmdxVTL5Rfs58aJuewrLOL5kTk0yFSLzVhKT0n+VivTNWplKtWIvirWYu7Oj6fNZ9G6HTx6TS+OPqJu2JFqpOJWpr8v0cp09Va1MpXEpwJRi41+ZxmvfPYNd597LKcfe0TYcWq8K/q0YuL1/Vi/Yy+XPqZWppL4VCBqqdcXrOPht5YyoFdLRp7aPuw4tcaJHZow/ZYTqZMaaWX6plqZSgJTgaiFFq/bwY9emEuP1g353WXdMLOwI9Uqxa1MOx1Zn5vUylQSmApELbNld2QYjbrpKYwd3EdjBoWkab10nlcrU0lwKhC1SOGBIm55No8NO/cxdkgOzerXCTtSraZWppLoVCBqkV//cyGfLt/C/7u8Gz1bNww7jvC/rUyvGqtWppI4VCBqiSmfrmTypyu56bT2DOjVKuw4UkpxK9PlG9XKVBKHCkQt8Onyzdz3j885vVNT7vquhtFIVKcfewTT1MpUEogKRA339ZZ8bp6SR9vGmfzpml4kJ+mMpUSmVqaSSFQgarDd+/Zz46RcDhQ5Tw/tS/06Gt6hOlArU0kUKhA1VFGR86MX5rJ0/U4eG9Sbo5pkhR1JKqBuegrjhuYw6Pg2PP7el9z6/Bz2Fh4IO5bUMhqsr4b647++4I3P1/PLCztzyjFNw44jh6G4lWmbRpn87rXFfLN9L08NyaFRVlrY0aSW0BZEDfTK/G949F9fMDCnFcNPahd2HKkEM+Om0/7bynTAmI/UylSqjApEDbNgzXbumDaXPm2zuf/SrhpGo4ZQK1MJgwpEDbJx5z5GTsolOzONJ67rQ3qKhtGoSdTKVKqaCkQNUbC/iJun5LElv4CnhuTQtF562JEkDkq3Mh39jlqZSvyoQNQA7s4v/raA3JVb+f2VPejaskHYkSSOiluZDujVkt+/GWllWnhArUwl9nQWUw0w8eMV/CX3a35wxtFc2L1F2HGkCqSnJPPIwB60CVqZrt2+hzGD1MpUYktbENXch19s4v5XFnF252b88KyOYceRKlSylenMr9TKVGJPBaIaW7FpN997bjZHN63LH67qSZKG0aiVSrcynb96W9iRpIZQgaimdu4t5IZJuSQZPD00h7rp2ltYm5VsZTrwSbUyldhQgaiGDhQ5tz0/lxWbdjNmUB9aN8oMO5IkgNKtTMd/+FXYkaSaU4Gohn7/5hLeWbyBey/uQv8OjcOOIwmkZCvTX6uVqVSSCkQ18/e5a3j8vS8ZdHwbBp/QNuw4koCKW5neeMpRPPPxCrUylcOmAlGNzPt6G3e9OJ/jj2rEvRd1CTuOJLCkJONnF3T+VivTDWplKhWkAlFNbNixl5GTc2laL50xg3qTlqKXTspXspXppWplKhWkT5lqYG/hAUZOzmPn3v08NSSHxnU1jIYcutOPPYIXblIrU6k4FYgE5+78dPpnzP16G48M7MlxzeuHHUmqoa4tv93K9Hm1MpVDoAKR4J7+4Cumz1nDj87uyLldjww7jlRjzRtkMG1Uf04+ugl3q5WpHIK4FQgz62Rmc0tcdpjZ7WZ2pZl9bmZFZpZTzjqSzWyOmb0cr5yJ7N0lG/jda4u4oFtzfnDG0WHHkRqgXp1Uxg3N4Vq1MpVDELef37r7EqAnRD7ogTXAS0AmcBnw5CGs5jZgEVDr9qss27CLW5+bw7FH1uf/ruyuxj8SMynJSTxwaVfaqpWplKOqdjGdCXzp7ivdfVFQPA7KzFoBFwBPxz1dgtmeX8jISbmkpSTx1NAcMtM0jIbEVulWppeN+YivNu0OO5YkmKoqEFcDUyu4zB+Bu4CDDnRvZiPNLNfMcjdurP5nZ+w/UMQPnp/D11vzeWJwH1o2zAg7ktRg53drznM3nsDOvfsZMOYjtTKVb4l7gTCzNOBiYFoFlrkQ2ODueeXN6+5j3T3H3XOaNm1aiaSJ4cHXFvP+0o3cf0lX+rZrFHYcqQX6tFUrU4muKrYgzgNmu/v6CixzEnCxma0AngfOMLMp8QiXSKblfs3TH37FsBPbcXW/NmHHkVpErUwlmqooENdQwd1L7n6Pu7dy93ZEdk+94+7XxSNcoshbuZWfvbSAk45uzM8vOC7sOFILlW5l+pO/qpVpbRfXAmFmWcDZwPQS1w0ws9VAf+AVM3sjuL6Fmb0azzyJ6pvte7hpch7NG9Zh9DW9SUnWz1MkHMWtTG898xheyF3NsAkz2b6nMOxYEhKrSZuROTk5npubG3aMCtlTcICBT37CV5t2M/2WE+nYrF7YkUQAeDFvNfdMn0+7xllMGN6XVtnqO1ITmVmeu0f9TZq+qobI3bnrr/NZsHY7f7yqp4qDJJTiVqbr1Mq01lKBCNGY977kn/PW8uPvduKszs3CjiPyP07s0ISXglamVz35qVqZ1jIqECF5e+F6fv/mEi7u0YKbT+sQdhyRMhW3Mu14ZD1umpLHC7O+DjuSVJHDLhDBAWg5DEvX7+S25+fQtUUDHrpCw2hI4ituZXrKMU25e/p8Xl/wTdiRpAqUWyDMrKWZ5QQ/eMPMjjCz3wJfxD1dDbR1dwE3TMwlMz2FsUP6UCc1OexIIockIy2ZJ67rTc/WDbl16lw+XrYp7EgSZwctEGZ2OzAX+DPwqZndQGTwvAygT7zD1TSFB4r43nOzWbd9L08O7kPzBhpGQ6qXzLQUxg/ry1FNsrhxUq4OXNdw5W1BjAQ6uXt/4FJgNHCOu//Q3bWNWUEPvLKIj7/czG8v60bvNtlhxxE5LMU/qMvOSmPYhFks27Ar7EgSJ+UViL3uvgXA3VcBSw5lfCT5X1NnruKZj1dw4ylHcUWfVmHHEamUZvXrMGXE8SQZDBk3g7Xb9oQdSeKgvALRysweLb4AzUv9LYdg5ldb+OXfF3Bax6bcfZ6G0ZCaoV2TLCZe34+de/czeNwMtuwuCDuSxFh5BeLHQF6JS+m/pRyrt+Zz85Q8Wmdn8ug1vUhO0hlLUnN0adGAp4fmsHrrHoZPmMmuffvDjiQxdNBONO4+saqC1ET5Bfu5cVIeBQeKeGpoDg0yUsOOJBJzx7dvzGPX9uamKXncNDmX8cP6kp6is/NqgoMWCDObAJQ1WJO7+4jYR6oZioqcO16Yx5J1O5gwvB8dmtYNO5JI3JzVuRkPXd6dO6bN4/bn5zL62t7aWq4Byutl+XKU61oDPwT0FeEg/vzOMl5bsI6fX3Acp3Ws/o2MRMpzeZ9WbNtTyP0vL+Tnf/uM3w7oph+BVnPl7WL6a/G0mbUHfgqcCjwIjItvtOrr9QXf8Ie3l3J571aMOPmosOOIVJkRJx/F1t0FjH53GdmZadx17rFhR5JKKG8LAjM7Fvg50Av4P2CUu+tIVBkWrt3BD/8yj15tGvLAgK76BiW1zh3ndGRLfgFj3vuS7Mw0bjy1fdiR5DCVdwxiGpFfTD9MZLfSAaB+8Yde8W8kJGLzrn3cOCmXBhmpPHmdhtGQ2snMuP+SrmzPL+SBVxfRMDOVK3Nahx1LDkN5WxB9iRykvhO4Ayj5ddgBfTUIFOwv4uZnZ7Np1z6mjerPEfXrhB1JJDTJScYjV/Vgx95C7p7+GQ0yUjmny5Fhx5IKOujvINy9nbsfFVzal5g+yt1VHALuzr3/+JyZX23hoSu6071Vw7AjiYQuPSWZJ67rQ9eWDfj+1Dl8unxz2JGkgsobrK+tmTUo8ffpZvYnM/th8eiuAlM+XcnUmau45TsduKRny7DjiCSMrPQUnhnWlzaNMrlhYi4L1mwPO5JUQHm/pH4ByAIws57ANGAV0BMYE89g1cXHyzZx3z8XctZxR3DnOZ3CjiOScLKz0pg8oh8NMlIZOn4myzdqcL/qorwCkeHua4Pp64Dx7v4wMBzoF9dk1cCqzfnc8txs2jfJ4g9X9SRJPwwSiap5gwwmj4h8ZAweN5N12/eGnEgORXkFouQn3hnAvwDcvShuiaqJXfv2c8OkWbjD00NzqFdHw2iIHEz7pnV5Zng/tu8pZPC4GWzL1+B+ia68AvGOmb1gZn8CsoF3AMysOVBrX92iIuf25+fy5cbdjBnUm7aN1X1V5FB0a9WAp4bksHJLPsMmzGK3BvdLaOUViNuB6cAK4GR3LwyuPxL4WfxiJbZH3lrK24vW88sLO3PS0U3CjiNSrfTv0Jg/X9OL+au3MWpKHgX7a/0OiYRV3mmu7u7Pu/sf3H1NiZvmAbXyk/Gf89Yy+t1lXNOvNUP6tw07jki19N0uR/Lg5d354ItN/OiFuRwoKmtMUAlTeae51jeze8xstJmdYxE/AJYDA6smYuL4bPV2fvziPPq2y+ZXF2sYDZHKGJjTmp+efywvz/+GX/59Ae4qEommvF9STwa2Ap8ANxAZrM+AS919bnyjJZYNO/cycnIujbPSefy6PqSllLd3TkTKM/LUDmzZXcgT//6Sxllp/EiniieU8gpEe3fvBmBmTwPfAG3cvVado7Zv/wFGTc5jW34hL97cnyZ108OOJFJj/OTcTmzdXcCj7yyjYWYa12sE5IRRXoEoPiiNux8ws9W1rTi4Oz97aQGzV21jzKDedGnRoPyFROSQmRkPDOjK9j2F/PrlhWRnpTKgV6uwYwnln8XUw8x2BJedQPfiaTPbURUBwzb+oxW8mLeaW888hvO7NQ87jkiNlJKcxB+v7smJHRpz57T5vLN4fdiRhPLPYkp29/rBpZ67p5SYrl9VIcPy/tKNPPDKQr7bpRm3n3lM2HFEarQ6qcmMHZJDlxb1uXnKbGZ+pW4CYdOR1jIs37iL7z83m47N6vHIQA2jIVIV6qanMGFYX1pmZzBi4iwWrq0VOyoSlgpEFDv2FnLDpFxSkpN4akgOWenlNt4TkRhpXDedySOOp256CkPGz2Tl5t1hR6q1VCBKOVDk3Dp1Dqs25zNmUG9aN8oMO5JIrdOyYWRwvwNFRVw3bgYbdtSqc2MShgpEKQ+9vpj3lmzkV5d04YT2jcOOI1JrHX1EPZ4Z3o/NuwoYMn4m2/MLy19IYkoFooTps1fz5PvLGXxCWwYdr2E0RMLWo3VDxg7OYfnG3Vw/cRZ7Cg6EHalWUYEIzFm1lbunf8YJ7Rvxy4s6hx1HRAInH9OEP13dkzmrtnLzs3kUHtDgflUlbgXCzDqZ2dwSlx1mdruZXWlmn5tZkZnllLFsHTObaWbzgnl/Fa+cAOu27+WmyXk0q5/OmEF9SE1W3RRJJOd1a84DA7rx3pKN3DltHkUa3K9KxO30HHdfQqQ1KWaWDKwBXgIygcuAJw+y+D7gDHffZWapwIdm9pq7fxrrnHsLD3DT5Fx279vP5BEn0ShLrbZFEtE1/dqwNb+Ah15fQsOMVO67uIsGzIyzqjp/80zgS3dfWXzFwV5YjwzrWNy4NjW4xOUrgzu0bpTJ904/mk5H1ovHXYhIjNx8Wge27i7gqQ++olFWOredpR+wxlNVFYirgakVWSDY6sgDjgYec/cZZcw3EhgJ0KZNmwoHy0hLZvS1vSu8nIhUPTPjp+cfx9b8Qv7w9lKys1IZ0r9d2LFqrLjvbDezNOBiYFpFlnP3A+7eE2gF9DOzrmXMN9bdc9w9p2nTppXOKyKJzcx48LJunHVcM+79x+f8fe6a8heSw1IVR2PPA2a7+2GNvuXu24B3gXNjGUpEqq+U5CRGX9uLfu0acccL83hvyYawI9VIVVEgrqHiu5eamlnDYDoDOBtYHPtoIlJd1UlN5qmhOXQ6sh6jpuSRt1KD+8VaXAuEmWUR+XCfXuK6AWa2GugPvGJmbwTXtzCzV4PZmgPvmtl8YBbwlru/HM+sIlL91K+TysTr+9G8QQbDJ8xi8ToN7hdLVpP6wObk5Hhubm7YMUSkin29JZ8rnvgYd/jrzSdqDLUKMLM8d4/6mzT9IkxEqr3WjTKZPOJ49u2PDO63cee+sCPVCCoQIlIjdGxWjwnD+7Jhx77I4H57NLhfZalAiEiN0btNNk8M7sOyDTu5cWIuews1uF9lqECISI1yWsemPDKwJ7NWbuH7z83W4H6VoAIhIjXORT1a8OtLuvL2og385MX5GtzvMKmXpojUSINPaMu23QU8/NZSGmam8YsLj9PgfhWkAiEiNdb3zziaLfkFjP/oKxrXTeN7px8ddqRqRQVCRGosM+MXF3RmW34h//fGEhpmpqpbZAWoQIhIjZaUZDx0RXe27ynk539bQMOMNC7o3jzsWNWCDlKLSI2XmpzEY9f2JqdtNrf/ZQ4ffLEx7EjVggqEiNQKGWnJPD20Lx2a1uWmyXnMWbU17EgJTwVCRGqNBhmpTBrRj6b10hn+zCy+WL8z7EgJTQVCRGqVI+rVYfL1x5OanMTgcTNZvTU/7EgJSwVCRGqdNo0zmXR9P/IL9jNk3Ew27dLgftGoQIhIrXRc8/qMH9aXtdv3MGzCTHbu1eB+palAiEitldOuEY8P6sPib3Zy4yQN7leaCoSI1GqnH3sEDw/swafLt3Dr1Dns1+B+/6ECISK13iU9W3LfRZ15c+F67pn+GTWp02Zl6JfUIiLAsJOOYmt+IX/61xc0ykrjnvOPCztS6FQgREQCt591DFvzC3jy/eVkZ6Ux6rQOYUcKlQqEiEjAzLjvoi5syy/kwdcW0zAjlav7tQk7VmhUIERESkhKMn5/ZQ+27ynkpy99RsPMVM7tWjsH99NBahGRUtJSknj8ut70bN2QW6fO5eNlm8KOFAoVCBGRKDLTUhg/rC9HNcnixkm5zF+9LexIVU4FQkSkDA0z05g0oh/ZWWkMmzCLZRt2hR2pSqlAiIgcRLP6dZgy4niSDIaMm8HabXvCjlRlVCBERMrRrkkWE6/vx869+xk8bgZbdheEHalKqECIiByCLi0a8PTQHFZv3cPwCTPZtW9/2JHiTgVCROQQHd++MY9d25sFa3dw0+Rc9u2v2YP7qUCIiFTAWZ2b8dDl3flo2WZuf34uB4pq7rhNKhAiIhV0eZ9W/OLCzry2YB0//1vNHdxPv6QWETkMI04+iq27Cxj97jKyM9O469xjw44UcyoQIiKH6Y5zOrIlv4Ax731JdmYaN57aPuxIMaUCISJymMyM+y/pyvb8Qh54dRENM1O5Mqd12LFiRgVCRKQSkpOMR67qwY69hdw9/TMaZKRyTpcjw44VEzpILSJSSekpyTxxXR+6tmzA96fO4dPlm8OOFBMqECIiMZCVnsIzw/rSplEmN0zMZcGa7WFHqrS4FQgz62Rmc0tcdpjZ7WZ2pZl9bmZFZpZTxrKtzexdM1sYzHtbvHKKiMRKdlYak0f0o0FGKkPHz2T5xuo9uF/cCoS7L3H3nu7eE+gD5AMvAQuAy4D3D7L4fuAOd+8MnAB8z8w6xyuriEisNG+QweQR/QAYPG4m67bvDTnR4auqXUxnAl+6+0p3X+TuSw42s7t/4+6zg+mdwCKgZRXkFBGptPZN6/LM8H5s31PI4HEz2JZfPQf3q6oCcTUw9XAWNLN2QC9gRhm3jzSzXDPL3bhx4+EnFBGJoW6tGvDUkBxWbsln2IRZ7K6Gg/vFvUCYWRpwMTDtMJatC/wVuN3dd0Sbx93HunuOu+c0bdq0cmFFRGKof4fG/PmaXsxfvY1RU/Io2F8UdqQKqYotiPOA2e6+viILmVkqkeLwrLtPj0syEZE4+26XI3nw8u588MUmfvRC9Rrcryp+KHcNFdy9ZGYGjAMWufsjcUklIlJFBua0Zlt+Ab99dTENMlL5zaVdiXzMJba4bkGYWRZwNjC9xHUDzGw10B94xczeCK5vYWavBrOdBAwGzihxmuz58cwqIhJPI0/twKjTOvDsjFX84a2lYcc5JHHdgnD33UDjUte9ROR019LzrgXOD6Y/BBK/vIqIVMBPzu3E1t0FPPrOMhpmpnH9yUeFHemgNBaTiEgVMTMeGNCV7XsK+fXLC8nOSmVAr1ZhxyqThtoQEalCKclJ/PHqnpzYoTF3TpvPO4srdP5OlVKBEBGpYnVSkxk7JIcuLepz85TZzPxqS9iRolKBEBEJQd30FCYM60vL7AxGTJzFwrVRf+oVKhUIEZGQNK6bzuQRx1M3PYUh42eycvPusCN9iwqEiEiIWjaMDO53oKiI68bNYMOOxBncTwVCRCRkRx9Rj2eG92PzrgKGjJ/J9vzCsCMBKhAiIgmhR+uGjB2cw/KNu7l+4iz2FBwIO5IKhIhIojj5mCb86eqezFm1lZufzaPwQLiD+6lAiIgkkPO6NeeBAd14b8lG7pw2j6IQB/fTL6lFRBLMNf3asDW/gIdeX0LDjFTuu7hLKIP7qUCIiCSgm0/rwNbdBTz1wVc0ykrntrOOqfIMKhAiIgnIzPjp+cexNb+QP7y9lOysVIb0b1elGVQgREQSlJnx4GXd2JZfyL3/+JwGGalc0rNlld2/DlKLiCSwlOQkRl/bi37tGnHHC/N4b8mGKrtvFQgRkQRXJzWZp4bm0OnIeoyakkfeyqoZ3E8FQkSkGqhfJ5WJ1/ejeYMMhk+YxeJ18R/cTwVCRKSaaFI3nUnX9yMjLZkh42by9Zb8uN6fCoSISDXSulEmk0ccz779kcH9Nu7cF7f7UoEQEalmOjarx4ThfdmwY19kcL898RncTwVCRKQa6t0mmycG92HZhp3cODGXvYWxH9xPBUJEpJo6rWNTHhnYk6OaZJGSFPuhOPRDORGRauyiHi24qEeLuKxbWxAiIhKVCoSIiESlAiEiIlGpQIiISFQqECIiEpUKhIiIRKUCISIiUalAiIhIVObuYWeIGTPbCKwMO0egCbAp7BDlSPSMiZ4PEj9joucDZYyFyuRr6+5No91QowpEIjGzXHfPCTvHwSR6xkTPB4mfMdHzgTLGQrzyaReTiIhEpQIhIiJRqUDEz9iwAxyCRM+Y6Pkg8TMmej5QxliISz4dgxARkai0BSEiIlGpQIiISFQqEIfBzM41syVmtszM7o5yexsze9fM5pjZfDM7v8Rt3c3sEzP73Mw+M7M6iZLPzFLNbGKQa5GZ3RPrbBXI2NbM/hXke8/MWpW4baiZfRFchiZSPjPrWeL1nW9mV8UjX2Uylri9vpmtNrPRiZYv+D/6ZvD/cKGZtUvAjA8Fr/MiM3vUzGLe0s3MxpvZBjNbUMbtFtz3siBj7xK3Vf594u66VOACJANfAu2BNGAe0LnUPGOBm4PpzsCKYDoFmA/0CP5uDCQnUL5rgeeD6UxgBdAupOdwGjA0mD4DmBxMNwKWB/9mB9PZCZSvI3BMMN0C+AZomEjPYYnb/wQ8B4xOtHzAe8DZwXRdIDORMgInAh8F60gGPgG+E4eMpwK9gQVl3H4+8BpgwAnAjOD6mLxPtAVRcf2AZe6+3N0LgOeBS0rN40D9YLoBsDaYPgeY7+7zANx9s7vHutN4ZfI5kGVmKUAGUADsiHG+Q83YGXgnmH63xO3fBd5y9y3uvhV4Czg3UfK5+1J3/yKYXgtsAKL+SjWsjABm1gdoBrwZh2yVymdmnYEUd38LwN13uXt+ImUk8l6pQ6SwpAOpwPpYB3T394EtB5nlEmCSR3wKNDSz5sTofaICUXEtga9L/L06uK6k+4DrzGw18Crwg+D6joCb2RtmNtvM7kqwfC8Cu4l8610F/N7dD/afM54Z5wGXBdMDgHpm1vgQlw0z33+YWT8iHyBfxjhfpTKaWRLwMHBnHHJVOh+R98k2M5se7Ab9PzNLTqSM7v4JkYLxTXB5w90XxSFjecp6DDF5n6hAxMc1wDPu3orIJuDk4E2ZApwMDAr+HWBmZyZQvn7AASK7Ro4C7jCz9iHkg8iH12lmNgc4DVgTZEsUB80XfIubDAx396JwIpaZ8RbgVXdfHVKuYmXlSwFOCW7vS2QX0LBEymhmRwPHAa2IfPCeYWanhJQxblLCDlANrQFal/i7VXBdSSMINufc/ROLHIhuQqSKv+/umwDM7FUi+xf/lSD5rgVed/dCYIOZfQTkENl/GUvlZgx2z1wGYGZ1gcvdfZuZrQG+U2rZ9xIlX/B3feAV4GfBZn88VOY57A+cYma3ENm/n2Zmu9z9fw7ShpRvNTDX3ZcHt/2NyP71cTHMV9mMNwKfuvuu4LbXgP7ABzHOWJ6yHkNs3iexPqhS0y9EiupyIt+wiw9sdSk1z2vAsGD6OCL7+I3IwaLZRA4ApwBvAxckUL6fABOC67OAhUD3kJ7DJkBSMP0A8OtguhHwVfBcZgfTjRIoXxqRgn97Avw/jJqx1DzDiM9B6so8h8nB/E2DvycA30uwjFcF798UIscf/gVcFKfXuh1lH6S+gG8fpJ4ZXB+T90nc/gPX5AuR3TJLiexb/llw3a+Bi4PpzkTOcJgHzAXOKbHsdcDnwALgoUTKR+Tb5LQg30LgxyE+h1cAXwTzPA2kl1j2emBZcBmeSPmC17cweF6LLz0TKWOpdQwjDgUiBq/x2UTO+PsMeAZIS6SMRIrYk8Ci4L3ySJzyTSVyjKOQyB6IEcAoYFRwuwGPBfk/A3Ji+T7RUBsiIhKVDlKLiEhUKhAiIhKVCoSIiESlAiEiIlGpQIiISFQqECJxYGbtyhqBU6S6UIEQiQ9D7y+p5vQfWCRGgq2GJWY2icgPITPM7KmgZ8CbZpYRzNfTzD4Nxu9/ycyyw00uEp0KhEhsHQOMAboQGSPnMXfvAmwDLg/mmQT8xN27E/n1670h5BQplwqESGyt9P8O0PeVu88NpvOAdmbWgEgDoX8H108k0hRGJOGoQIjE1u4S0/tKTBcPYy1SbahAiFQhd98ObC3RO2Aw8O+DLCISGn2jEal6Q4EnzCyTyHDTw0POIxKVRnMVEZGotItJRESiUoEQEZGoVCBERCQqFQgREYlKBUJERKJSgRARkahUIEREJKr/Dy9Zkc6Pv9yPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.array(all_results[\"rho\"][:]), np.array(all_results[\"mseF,PSNR (p2plane)\"][:]))\n",
    "plt.title(\"optimal rho for D2 distortion\")\n",
    "_ = plt.xlabel('rho')\n",
    "_ = plt.ylabel('RSNR')"
   ]
  },
  {
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}