igl_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b9849bc6",
   "metadata": {},
   "source": [
    "## Installation "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6913257",
   "metadata": {},
   "source": [
    "activate environment and install \n",
    "<br>\n",
    "conda install -c conda-forge igl \n",
    "<br>\n",
    "conda install -c conda-forge meshplot\n",
    "<br> \n",
    "pip install numpy\n",
    "<br>\n",
    "pip install scipy\n",
    "\n",
    "<br>\n",
    "https://drive.google.com/drive/folders/1X34dBiyU0x6Pz-kbsuseJ6mRlnzZ2xF-?usp=sharing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6565aa1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import igl\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "from meshplot import plot, subplot, interact"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bada2196",
   "metadata": {},
   "source": [
    "### Loading a mesh \n",
    "\n",
    "In order to load a mesh igl comes with pre-build functions that can read and parse mesh files. <br>\n",
    "A .mesh file contains the list of vertices and triangles, vertices are the 3d coordinates of points and triangles stores the index corresponding to each vertex in a triangle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5eb446d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#open bunny.mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "15b3c8c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "v,t, f = igl.read_mesh('./bunny.mesh')\n",
    "# v stores the vertices and f stores the faces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ef719731",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5433, 3) (6966, 3)\n",
      "## Vertices ## [[-0.0260146   0.112578    0.0363871 ]\n",
      " [-0.0321783   0.174119   -0.00263321]\n",
      " [-0.080718    0.152855    0.0302446 ]]\n",
      "## Faces ## [[2784 2497 2027]\n",
      " [1077  225 1060]\n",
      " [ 425  450  381]]\n"
     ]
    }
   ],
   "source": [
    "print(v.shape, f.shape)\n",
    "print('## Vertices ##', v[:3])\n",
    "print('## Faces ##', f[:3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9ef3aa1",
   "metadata": {},
   "source": [
    "### Mesh statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6be3a2c3",
   "metadata": {},
   "source": [
    "![alt text](https://www.researchgate.net/profile/Marco-Attene/publication/225913168/figure/fig2/AS:302290532618241@1449083087214/Meshes-Irregular-semi-regular-and-regular.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feb0f490",
   "metadata": {},
   "source": [
    "Libigl contains various mesh statistics, including face angles, face areas and the detection of singular vertices, which are vertices with more or less than 6 neighbours in triangulations or 4 in quadrangulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a80d0ece",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total vertices :  5433\n",
      "Irregular vertices  no :  2139 39.37%\n",
      "Area average: 8.356721028933905e-06 Min:3.1288891820704787e-07 Max:9.39460837031244e-05 Sigma:0.5728664262127175\n"
     ]
    }
   ],
   "source": [
    "irregular = igl.is_irregular_vertex(v,f)\n",
    "v_count = v.shape[0]\n",
    "irregular_v_count = np.sum(irregular)\n",
    "irregular_ratio = irregular_v_count / v_count\n",
    "print('Total vertices : ', v_count)\n",
    "print('Irregular vertices  no :  {} {:.2f}%'.format(irregular_v_count, (irregular_ratio * 100.)))\n",
    "\n",
    "area = igl.doublearea(v,f) /2.0 # igl.doublearea computes twice the area for each input triangle\n",
    "area_avg = np.mean(area)\n",
    "area_min = np.min(area)\n",
    "area_max = np.max(area)\n",
    "area_ns = (area - area_avg) / area_avg\n",
    "area_sigma = np.sqrt(np.mean(np.square(area_ns)))\n",
    "print('Area average: {} Min:{} Max:{} Sigma:{}'.format(area_avg, area_min, area_max, area_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca06e6c9",
   "metadata": {},
   "source": [
    "The first row contains the number and percentage of irregular vertices, which is particularly important for quadrilateral meshes when they are used to define subdivision surfaces: every singular point will result in a point of the surface that is only C^1.\n",
    "The second row reports the area of the minimal element, maximal element and the standard deviation. These numbers are normalized by the mean area, so in the example above 5.33 max area means that the biggest face is 5 times larger than the average face. An ideal isotropic mesh would have both min and max area close to 1.\n",
    "\n",
    "The third row measures the face angles, which should be close to 60 degrees (90 for quads) in a perfectly regular triangulation. For FEM purposes, the closer the angles are to 60 degrees the more stable will the optimization be. In this case, it is clear that the mesh is of bad quality and it will probably result in artifacts if used for solving PDEs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "d820dd7f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(6966, 3) (6966, 3)\n"
     ]
    }
   ],
   "source": [
    "angles = igl.internal_angles(v,f)\n",
    "print(angles.shape, f.shape)\n",
    " #f by poly-size array of internal angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e71cbf90",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "f51b3cfc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d6eadaac7df5444a942da1d075131fe9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016840…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb1789f8fd0>"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(v, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "badaaf57",
   "metadata": {},
   "source": [
    "### Mesh Decimation\n",
    "\n",
    "Given a high resolution mesh with too many triangles, find a “well approximating” low resolution mesh with far fewer triangles. One family of mesh decimation methods operates by successively remove elements from the mesh. In particular, Hoppe advocates for successively remove or rather collapsing edges. \n",
    "![alt text](https://www.intel.com/content/dam/develop/external/us/en/images/23792-figure1-165263.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "38844a12",
   "metadata": {},
   "outputs": [],
   "source": [
    "v, f = igl.read_triangle_mesh('./bunny.obj')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "3dbb57ba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((34834, 3), (69451, 3))"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v.shape, f.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "4fd33c26",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(34834, 3) (69451, 3)\n",
      "(999, 3) (999, 3)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ef8424982254204b8c9afb77ec92a33",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016663…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb17877fe10>"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(v.shape,f.shape)\n",
    "_,v_,f_,_,_ = igl.decimate(v,f, 1000)\n",
    "angles = igl.internal_angles(v_,f_)\n",
    "print(angles.shape, f_.shape)\n",
    "plot(v_, f_, angles)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "288b2cae",
   "metadata": {},
   "source": [
    "### Mesh Subdivision\n",
    "Given a coarse mesh (aka cage) with vertices V and faces F, one can create a higher-resolution mesh with more vertices and faces by subdividing every face. That is, each coarse triangle in the input is replaced by many smaller triangle. <br>\n",
    "\n",
    "\n",
    "#### \"igl.upsample\"\n",
    "adds new vertices at the midpoint of the edges. Each resutling triangle thus is split into four triangles.\n",
    "#### \"igl.loop\"\n",
    "computes the weight combinations of vertex neighbours and moves the vertices of the mesh to those locations, while performing smoothing as it is refined. The combinatronics are the same as upsample however smoothing results in convergence to a C1 surface. \n",
    "\n",
    "#### \"in\" plane\n",
    "\n",
    "This method will not change the point set or carrier surface of the mesh, new vertices are added on the planes of existing triangles and vertices surviving from original mesh are not removed.Thus, the combinatronics of the mesh are changed. This is done using the \"subdivision\" rule."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95ccc17d",
   "metadata": {},
   "source": [
    "![barycenter](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAATIAAAClCAMAAADoDIG4AAABs1BMVEX///9vWX+vprcAAAD/Vx7//P//+///OgD/+P//QQD/NwD/RQD/WRzhwvT7+vz/QgC+h98sLCzIyMje3t6QkJD/VwDr2fn/0dL98///9PoAAD8AADn/rrD/5u45AF20tLROTk7/1t3btfESEhLU1NRZWVnx8fHz6v8XAEYAADtzRo9KAG3/Vy3/7e//wcf25v//ioEvAF8+AGfRr+TcyOnl1u//XTYmAE//m5anhbv/r6n/ubv/d2eKV6b/29//jHr/ZEuHYaGdacDXw+PJrdhnZ2erq6tBQUEdHR19fX2goKAtLS3kz/ujAGLfMQDSJwyjRJSrJn7RJByxAIbDC03oPB6/LkifYIK0VmvedIWNVbNsK4//bmnMtNyukMKVZ7F2Q5JhMH1UFXL/dFhYGXpkLYOqcsexg82jX82/mNt/QqXGo92tb9GAVpmqib29os99SZvQMjq0TYGlTZC1fbi8AEGvbbu/HTawI2f/lIuXcq20jc3vvNHznaWfLGrQp+i4PW2wG3WBNKjNy+q4aZuvYq69AC7fd5azPVXXMDJnQ35rX4OyIFWIG2mYLInPFRe3WrkYBkS9AAALcUlEQVR4nO2d/UPaSBrH5zYJb8FqVQQUtNQXqCAg8lZAZKtWbal9uevdbq26ispRa+1u7Yt6bfdeerfddvfuT74EBCGZZAZrTQbm8wMmYUKGr09m5nnmyQAAhUKhUCgUCoXSfXNE6yqQxmXmitZVIIwRZpEZ07oSZLHEdDO3tK4EWTCjYPi21pUgim5mqWee6dG6GiQxz4j0al0Nghgrd5fDl7WuB0F03x4SXntv0zuTQtEbV/Na14Aw+u4s373n0LoWJGG6b7Na+QOtq0ESYw+sLGv7o9bVIAnHn0TJ7mldDaI4eGiz/nkypnU1iML/l+++v5GgmuGzmvB64uCRj2qGzfge4Ixu8GiyT+uakELsWhSATApQzbBZXxFeXI8B1QwXR2JVeM2aDcLreIRqhsGaXxSLC9jFnT1qZ2hM/rXy33Cy/Gc8EtKyNkSw6qs45MFcZX+PaoZia73y12GcqGyM+6lmquT7qwPYXPBkg96b6uz9UN1Khk82TBvUzlTw+marmx5zdcu0Qe1Mmev+2qaTd1c3HRubNEirADe5e7pTcNU2TRvbVDM4u4N1yrgKp9umjQ2qGZSto7odN+s83TFtUc1grPoamnmzp27Hu0XbMwh74w27Jz7TCV4/tTMZsf7GGd+qz3SC17/NXWBtiGB9u3F/wtJoViH/ONWsAUdiTnIkl27cF+zswmpDBGsRqQ3VfKYqoUjZzmI0iFaG8z+THrIHpCKGIhsARHfo/VmmGiirw9SZlR4S2rOrReAoXkyddM7KuvzYY5fskDeyArx/NVxAhXRPHjbTm8rIDu2X/E92LqA+OscrjMe+hXWGbqO01TqeAlOT0xdRKX3D7e68T0RhbwQ8kiNeMFX0Jp5cRK30jGhJR5EuL+StmaT0SL4kDDISj75+rXSN42kxNLkPijvy4VYwDj2DagZM64ldAAxyt9t9GpptIJaA9K/tg6EYAxtHfU9LsEBFPAg5CETN2roPcLy/IQbKvLB8suSMwkkxX1vbGdgb34HaGACeuJJzFO1vZ836+qMgVISNMoDBCG/MBGbb2c6kgbJ6xNw8BaK+61+hMkTgGJQGyuoo5+Yp0L6a7U+qBHOyvEn5zVmfLGDUHvjVbIVj7Srvzg60pZ2tJmCuUo2wzGeqZ3egHe1s5VvVtyXzTFKm2tDO+nzQ0UUNh2VC9f2pa/vnWR0SgAbK6skp+ExVpvrbTDP5VJwU2TyTlHbT7JkfNV1kNyMKgFJbaWaIrKGKOFlFn6lKW2k2pz7CKFOQzzNJaSfNtjFc6/rcPCVKvl10oZYgeg0jcdjNO9GFSr6pL68OCUw/xynVKZ1nglHqbws7C53mrKuh7jNVud7f+nZmKE1vYhVE+ExV1kQ7Q/cmBOO4v8z/eIxTcsKIl/L5zHf80+cXLSzaAW+12l5iiSHNzVPi2QOb1XandVNqf7KxrPUu1oJISJ/phAPxI3ms9pFIBCtjra8fBzHy6+wBvOyoV2XJVr+wYvrF+26Z/3icCrAppEPEGWW5eRCcyYevhRvzTevemLFEcV+8LYOFzjBq4AXJzZPCJfmCe/Xd698nO86lenrkaKu65T7kc2nVET4kN68RLs3nypMEDnCEiFmSi8NXFyhzpgN8UuX+lOfmNRIMBE4DkUeJFs3YPvY3zrbZM8aMXUkYQ1zt1rXH2YZRyBE0v494DH5ZoGwiaY6nFGYt5bl5NTwFi0vSoe75WtHOVhOwji1YsMxAO8dgQOFz3BlzUj4CmW7F9kwpUJY9tMDmR9w8dJ7JPWM8hDaBz1tvKa9ov+JXcrpgXQEsN895aAwr9RnTg63ma07/oPZuuStoPCTPzTOlzBmVEe7zFtMslEA4gu4kH3DVD9U8ksaMc/E59QFwi9lZ3eIOSpjSOePhqRUZzA23YJpVyqI9ZaWVNDMNYgWds+HOQk2Y+tw8e47FCAdxraTZ/qRK0lg9zlQg4Kr0la6az2TP8Sm8yEYLaeY/Qpc5gQsWzGWvvZqbl80YkxgzTpWTV/wtEtaQLO6A4sRrZ/8mmIw7bDzEFUzkiV84KUa+bntYU3F1iF7723cvPx8kjRlkcK2RFb8pSvbDiKb3DhA7izfjeWOzWpd/blIw4YLjkSKIET1dFysWj86yDsGlH8XFxv/e/FPlXRtbHTuYvY0uEdwkb2K6o/kHnYMfBcmsL/kwTki7jtIq2IwQ3W9OPZ07HhRal/fNnRaMx1/wVuvdubLX3oTgxQ4QLW5skd3+5yP/8FYeW8XGk+NdwFT8+aU4AIZ77QpwITBXBKZtgjUrToFZX76000z776kOxAzV3Dx7xpLByWwRmSoJL6ZNgjWLdomroHL447KJGXNtIFbOzbs0tHSrZ0Lw2tNYbXqlkJdkzUK+riL+dPZEQ0QsVQCgh2GuLDK9Uq8dBcmarW8DrgMzE8yU7GyIiIm5eVeYEfGn+sTfT8uGzQVkOKOKNzKOLqRLvAns6X/OZSxIWqxOzxizJG70VH6mVewK0HPtJ1feHCdzDZe1CG7JNBuQPfA1kxxhhsCl3t7e6q/0GYI5I3KuvYJ3kEg74yYxM6eDgQAkIhbM9TC3wKXLC8zQ6UHBa4+rz7WfIGhGoBcwhxfAEgdisJHbhNHLLADxl1qH6g87XXGsoVofiXa2hRMoyxZ4l4I5xO29zOjQLaZRMiAO3szSCRYIBGoWxQiUuWcsyiHEZBgIejFDV6SSlSdY4i7U/dk3qP4ko/7YQ65v5AybwyoPE6rm5nHBXCfKaydNsxAqUOZMGtW/M2dRlyQbNiK8dsI0QwTKDCleOhCTgczNc7pYNqn20Gt+kiDNHINqIwwuzebQDTgyNw8gvXaSNNsdVBkWNSTVKZO14FzJnTTHlb32PDnLUm0qPxdhj7PobFgRQxxtieVyoteuNFTLJ25gfYjmKE/FZTPGFO5Mm0pungRP2KLktZOi2bZCE+JuYi5Xcd08KM4U25mCdgVkaKawuMNE2KKYIwZDITdPAUOwAPfaZ/sJ0Ay6fJQzaW52Lhed8dNI2WuXD9VmfbrXDLa4gyHViRyIyVBcN08RuNc+69P7qtHyQFnt6YbmkObmYWHP8DKvXfeaTUpz1oMs1kBMhiQ3DxeI1z7rw3oaVCumJIs72HN8+oxL0Kusm6eK4LVL5tpXB/RsZ9sNgbJsATepDoILw2dSQOq1z+lYs9hA3eMe4tMNzeSIScjyXzC9VvbaT8+fG9DtvVkXKHMfWg6bGVrJ4c/QadTR4LXPDeh0GZdQf3W219R8Up0MvDUgVHAnLbWhml41O9qqNPXnIRj2GhBqcOl41WvXp2aOSKVWWAMxA7JfmECsm3cCYgKu5rUf61Gz/YRoZEGMpxsEOq4ii6DWzavQdQlRwPmW5cUMeR1qlt+aFu8mzIEYhmRvv/8nxhAFKVnFa5/xgJJv518f9PRk3auP1o8HBfhcLgS0ZPk37PJ99BQyhmRAzIoXWosnyzbbv/XzBGfHXStrffAf/PJIyd7ZrFbbB+QHdWFe0JDOlVOXX2CW/2rMj46OzncLG9/wLGtl17/B5foRosBaOdH4F+QHTWNf8pn4ibbvhMr2DA2NaSbZ8OL8/AIjaLa7LFjZwz9g8+uviAKfPotf8LdPqA/awL7kJ9Fu+R0AlhiBbs0kmwdgjOkV7P533tZ5D/9E9I25+1+b7SV6rSDcG1Mgf+fhg3scuMnMj/UsMlrZ2eLCzZvzYsYhMOz+b7+JoAVGjxn98ApjhQK85r+CoUP8F1xmhJfupR5U6a/Eomjjw2e4OoZkeDQjWQXmyjld+myIN+aZ6tDuki00f+L5SdZ0RG5BvDFvLoycUwWaZZhZWJBnzumbIWZ0ZIgZ1uryvUtLS71a/b/OSq/Q/i5q1foTylg3af9lCoVCoVAoFAqFQqFQKBQKhdIS/B+0yIR/FtDkrgAAAABJRU5ErkJggg==)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "821efc86",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46646cffe782411eafa61ad5eb5da625",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016840…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cf499a9649464a29b8c78dbf2abb54ac",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016663…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "94dd3e365c76440f870e606069c61c14",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016973…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "39590c155f6441d1a16727cf831c4caf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016663…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb178a18be0>"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uv, uf = igl.upsample(v_, f_)\n",
    "lv, lf = igl.loop(v_, f_)\n",
    "fv, ff = igl.false_barycentric_subdivision(v_, f_)\n",
    "ia_l = igl.internal_angles(lv, lf)\n",
    "ia_u = igl.internal_angles(uv, uf)\n",
    "ia_f = igl.internal_angles(fv, ff)\n",
    "plot(v, f,  shading={\"wireframe\": True})\n",
    "plot(uv, uf, ia_l )\n",
    "plot(lv, lf, ia_u )\n",
    "plot(fv, ff, ia_f )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd8b3a28",
   "metadata": {},
   "source": [
    "### Delaunay Triangulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c15d95ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 3 0]\n",
      " [3 1 0]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.spatial import Delaunay\n",
    "points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])\n",
    "                    # 0    #1          #2        #3\n",
    "tri = Delaunay(points)\n",
    "print(tri.simplices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "705e150b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 3, 0],\n",
       "       [3, 1, 0]], dtype=int32)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tri.simplices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aac2dc53",
   "metadata": {},
   "source": [
    "### in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7305718a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "78b2dcab",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3caa022c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXiU5b3/8fedfSEkrAlb2HcVgQCuVVyQTaytrYrVavVQW7Uce+qxx7YSsSoutaKoSCkuba2nP4+twYCIK1ZEQNzIAoSwJRDCmn2SzMz9+yMBhyxkCJNMZubzui4uMvM8zHwfknzuuZ/l+xhrLSIiEvjC/F2AiIj4hgJdRCRIKNBFRIKEAl1EJEgo0EVEgkSEv964e/fudsCAAf56exGRgPT5558ftNb2aGqZ3wJ9wIABbNy40V9vLyISkIwxu5pbpl0uIiJBQoEuIhIkFOgiIkFCgS4iEiQU6CIiQSKgAj0zP5Mpr0/hrJfPYsrrU8jMz/R3SSIiHYbfTls8VZn5maSvTcfhcgCwr2If6WvTAZgxaIYfKxMR6RgCJtAXblp4PMyPcbgcLFj3JAnOiXSOjSQhJoKEmAg6x0QSHRGGMcZP1YqItL+ACfSiiqImnz9SU8xNy9Y3ej4y3JAQ823IJ0RH0jk2wuO5SDofW9boubrHMZHhbb1ZIiI+EzCBnhKfwr6KfY2e7xGbzLKfnkuZo5Yyh5MyRy2lDufxr8scTkrr/95xsKL+eSfl1c4W3zMqPOzbAaE+5Dt7hL/njOCEdTxmC9ERGhREpH0ETKDPHTf3hH3oADHhMfxqwt1MHNj1lF/P5baUV3uEflX9gFBdezz0Sx0eX1fVUuao5UBZ9fHBoqLG1eL7REWEnfCpv7nZQt3A0PTAERURUMeuRcRPAibQjx34XLDuSY7UFNMjNplfTbi71QdEw8MMibGRJMZGtroml9tS7hH8pR6zhBNnCJ4zh1r2l1YfX1bpxaAQHRHW5C6ipmYLDXcbHZstRIZrUBAJdgET6FAf6uXjuOPVTbxx93cYlpzg13rCwwyJcZEkxrV+UHC63PUzhboBobSq8YBQVn3irqTSqlr2lVQdnz1U1bY8KMREhjV9/KCZ2YLngHHs7wgNCiIdWkAFejCKCA8jKS6KpLioVr9GrctNeYPdRM3NFo4PHA4nhUerji9z1LpbfJ/YyPBmjyl0jo0kIbrxQeYTZxIaFETakgI9CESGh9ElPoou8ac3KDQ8kHzibMFzxvDtrqTCI1XHdyVVO1seFOKiwpsJ/ca7i5qaLXSK1qAg0hwFugB1g0LX+Ci6nsagUON0NxoASluYLZRU1rDncOXxdWu8GBTio8IbHTfw3JXkeZZRQnTj2UKnmAjCw3SNggQfBbr4TFREGN06RdOtU3SrX6Pa6Wo8I6ifLZQ2M1s4UlnD7sOVx89UqnG1PCh0Or57qPkZQXOzhYSYut1LYRoUpINRoEuHEh0RTnSncLqfxqDgqHU1mgmcsCupidnCofIadh6sOL6s1mVbfJ/Gxwya2JUU2/QFbAkxEXSK0qAgvqVAl6ATExlOTGQ4PRJaNyhYa6l2upueEXgcW2h4AduB8mry6y9eK62qxek++aBgTN1MoXOTA0ITB56bmC3Ea1AQDwp0kQaMMccHhZ6tPDPWWouj1n3C9QcnHlto6gK2WvaXOsgr/nbAcJ3ioND04HDiqakNdyXFR4Wr71GQUKCLtAFjDLFR4cRGhdOzc+tew1pLlcfuo4YzAs/ZgucFbPtKHGwt/vYK6BbGBMKODQqxkc2EftMXs3muE6dBoUNQoIt0UMYY4qIiiIuKILlzTKtew1pLZU3DQaH52cKxdQqPOihzlB1/rqVBITzMHD/Q3Dj0vZstxEZqUDhdCnSRIGaMIT46gvjoCFISWz8oVNS4TpwZtDBbKHM4KThSeXxZebWzxUEhIszQqbmGd17OFmIiQ7ttdouBboxZBswEiq21ZzSx3AALgelAJXCztXaTrwsVEf8wpu7Td6foCHoltu413G5LRY2ThgeZG88OTnxcd41C3Xrl1U6sF4PCyfodNTx+cGLbi7a/l0JmfiYLNy2kqKKIlPgU5o6b69Mb9HjzCf0lYBHwSjPLpwFD6/9MAp6v/1tEBICwsGP3J2h936Njg0JTp502fM5ztrDrUKXHdQstt81ueC+FxmcXtdwor6l7KbTHXddaDHRr7RpjzICTrHIV8Iq11gLrjDFJxphe1trGzctFRFrpxEEhtlWv4XZbyms82mU3HAgaHEs4tuygD+6lsDXqcappfNe1hZsWtl+ge6EPsMfjcUH9c40C3RgzB5gDkJqa6oO3FhHxXliYqb8K+DTbZlc3cdygiXspfDtw1OKIPERTe3Kauxtba/gi0Jva2dTkni5r7RJgCUBaWlrLl+KJiHQwJ9xLoUvL6+cVl/PA8ixsbRIm6mij5SnxKT6rzRdt6wqAfh6P+wJ7ffC6IiIBq8xRy8Mrcpj61Bq+3HOUGX1/Qkz4iWcaxYTHMHfcXJ+9py8+oWcAdxpjXqPuYGiJ9p+LSKiy1vKvLwt5eEUuB8ur+eH4ftwzdTjdO0WTmd+Te99fgIksoZc/znIxxvwduBjobowpAOYBkfWFLwZWUHfKYh51py3e4rPqREQCyObCEtIzsti46whj+iWx9KY0xvRLOr58xqAZ3L0sjJ+cP5BfTxvh8/f35iyX61tYboE7fFaRiEiAOVJRwxPvbOHV9bvpGhfFY9ecxTXj+rZ74zRdKSoi0kout+XV9bv5wztbKHM4ufm8AfznZcNO6+bzp0OBLiLSCht2Hmbem1lk7yvl3EHdSJ81muEp/r1xvQJdROQUFJc6eGRlLv/8opDeiTE8O3sc089M6RA9ZBToIiJeqHG6efGTHTz93jZqXZY7Jw/h55MHExfVcWK041QiItJBfbT1AA8szyL/QAWXjezJ72aOon+3eH+X1YgCXUSkGXsOVzL/rWxWZ+9nYPd4XrxlApOH9/R3Wc1SoIuINFBV4+L5j7az+KPtRIQZ7p06gp9cMIDoiMZdFDsSBbqISD1rLW9vLuL3mTkUHq1i1pje3Dd9ZKtvDtLeFOgiIsC2/WWkL8/ik7xDjEhJ4H/nnMOkQd38XdYpUaCLSEgrddSy8N1tvLx2J3FR4cy/ajSzJ6YSEe6L3oXtS4EuIiHJ7ba88UUhC1bmcqiimusmpHLPFcPpGh/l79JaTYEuIiHnm4IS7s/YzBe7jzI2NYkXb57AmX1becPUDkSBLiIh43BFDY+vyuW1DXvoFh/NEz8Yw/fG9mn3JlptRYEuIkHP6XLzt8/qmmhV1ri49fyB/OKyoad1K7qOSIEuIkHts/xDzMvIIreojAuGdCd91iiG9PRvE622okAXkaBUVOLg4RU5ZHy1lz5JsTx/wzimntExmmi1FQW6iASVaqeLP/97B4vez8Pptvzi0qH87KLBxEZ17Ks8fUGBLiJB44PcYua/lc2OgxVMGZXM72aOol/XOH+X1W4U6CIS8HYdqmD+8mzeyy1mUPd4Xv7JRC4a1sPfZbU7BbqIBKzKGifPfbCdJWvyiQw3/M+0Edxy/kCiIgLvKk9fUKCLSMCx1pL5zT4ezsxhb4mDq8f24dfTRpDcOTCaaLUVBbqIBJQtRWWkZ2Txaf4hRvXqzNPXjyVtQFd/l9UhKNBFJCCUVNXy1LtbeeXTXSTERPD7757B9RNTCQ+Sqzx9QYEuIh2a2215/fMCHn07l8OVNcyemMqvpgynSwA30WorCnQR6bC+3HOUeRlZfLXnKGn9u/DyrImc0Sfwm2i1FQW6iHQ4B8ureeztXP6xsYAeCdH88doxfPfsPkF9lacvKNBFpMNwutz8Zd0unly9laoaF3O+M4i7LhlCQpA10WorXgW6MWYqsBAIB5Zaaxc0WJ4I/BVIrX/NJ6y1L/q4VhEJYp9uP0R6RhZb9pdx4dDuzLtyNEN6dvJ3WQGlxUA3xoQDzwKXAwXABmNMhrU222O1O4Bsa+2VxpgewBZjzN+stTVtUrWIBI29R6t4aEUOmV/vo2+XWF64cTxTRiVr90orePMJfSKQZ63NBzDGvAZcBXgGugUSTN13oBNwGHD6uFYRCSKOWhdLP87n2Q+247aWuy8bxk8vGkRMZPA30Wor3gR6H2CPx+MCYFKDdRYBGcBeIAG41lrrbvhCxpg5wByA1NTU1tQrIkHgvZz9zH8rm12HKpl2Rgq/mTGSvl1Cp4lWW/Em0Jua99gGj68AvgQuAQYDq40xH1trS0/4R9YuAZYApKWlNXwNEQlyOw5WMH95Fh9sOcCQnp34662TuGBod3+XFTS8CfQCoJ/H477UfRL3dAuwwFprgTxjzA5gBLDeJ1WKSECrqHay6IM8/vzxDqIiwvjtjJH8+LwBRIaHZhOttuJNoG8AhhpjBgKFwHXA7Abr7AYuBT42xiQDw4F8XxYqIoHHWsvyr+uaaBWVOvj+uL7cO204PRNCu4lWW2kx0K21TmPMncAq6k5bXGatzTLG3F6/fDHwIPCSMeYb6nbR3GutPdiGdYtIB5ezr5T0jCw+23GYM/sk8uwN4xjfv4u/ywpqXp2Hbq1dAaxo8Nxij6/3AlN8W5qIBKKSylqeXL2Fv6zbRWJsJA9ffSbXTuinJlrtQFeKiohPuNyWf2zcw+OrtnC0soYfndOfX14+jKQ4NdFqLwp0ETltm3YfYd6bWXxTWMLEAV1JnzWaUb07+7uskKNAF5FWO1BWzaNv5/L65wUkd45m4XVnM2tMb13l6ScKdBE5ZbUuNy+v3cnCd7fhcLq4/aLB3HXJEOKjFSn+pP99ETkln+QdJD0ji23F5Vw8vAf3zxzFoB5qotURKNBFxCsFRyp5KDOHlZuLSO0ax9Kb0rh0ZE/tXulAFOgiclKOWhdL1uTz3Id5APxqyjBuu1BNtDoiBbqINMlay+rs/TyYmc2ew1XMOLMX980YSZ+kWH+XJs1QoItII9sPlPPA8mzWbD3AsOROvHrbJM4boiZaHZ0CXUSOK6928sx721j2yQ5iIsK5f+Yobjy3v5poBQgFuohgreXNL/fy8Iocisuq+WFaX/576gi6d4r2d2lyChToIiEua28J6RlZbNh5hDF9E3nhxvGMTVUTrUCkQBcJUUcra3jinS28+tluusRF8ej3z+QH4/sRpiZaAUuBLhJiXG7L39fv5ol3tlDmcHLTuQO4+/JhJMZG+rs0OU0KdJEQsnHnYeZlZJG1t5RJA7vywFWjGZGiJlrBQoEuEgKKSx0sWJnLG18U0isxhmeuH8vMs3rpKs8go0AXCWI1Tjcvrd3B0+/lUeN0c8fkwdwxeQhxUfrVD0b6rooEqTVbD5C+PIv8AxVcOqInv5s5igHd4/1dlrQhBbpIkNlzuJLfZ2azKms/A7rFsezmNC4ZkezvsqQdKNBFgoSj1sXzH25n8UfbCTOGe64Yzm0XDiQ6Qk20QoUCXSTAWWtZlVXEg2/lUHi0iivH9Oa+6SPolagmWqFGgS4SwPKKy0jPyObfeQcZkZLAa3PO4ZxB3fxdlviJAl0kAJU5ann6vW28+MlO4qLCSb9yFD86pz8RaqIV0hToIgHE7bb884tCFrydy8Hyaq5N68c9Vwynm5poCQp0kYCxubCE+9/czKbdRzm7XxJLb0pjTL8kf5clHYgCXaSDO1xRw+OrtvDaht10i4/i8WvO4vvj+qqJljSiQBfpoJwuN6+u380f3tlKebWTW84byH9ePpTOMWqiJU3zKtCNMVOBhUA4sNRau6CJdS4GngIigYPW2ot8WKdISFm/o66JVs6+Us4b3I30WaMZlpzg77Kkg2sx0I0x4cCzwOVAAbDBGJNhrc32WCcJeA6Yaq3dbYzp2VYFiwSzohIHj6zM4c0v99I7MYbnbhjHtDNS1ERLvOLNJ/SJQJ61Nh/AGPMacBWQ7bHObOANa+1uAGttsa8LFQlmNU43yz7ZwdPvbcPptvzikiH87OIhxEbpKk/xnjeB3gfY4/G4AJjUYJ1hQKQx5kMgAVhorX2l4QsZY+YAcwBSU1NbU69I0PlwSzHzl2eTf7CCy0Ymc//MUaR2i/N3WRKAvAn0puZ6tonXGQ9cCsQCnxpj1llrt57wj6xdAiwBSEtLa/gaIiFl96FK5r+Vzbs5+xnUPZ6XbpnAxcO1t1Jaz5tALwD6eTzuC+xtYp2D1toKoMIYswYYA2xFRE5QVePiuQ/zeGFNPpFhhl9PG8FPzh9IVISu8pTT402gbwCGGmMGAoXAddTtM/f0JrDIGBMBRFG3S+aPvixUJNBZa1nxTREPZWazt8TBd8/uzf9MH0ly5xh/lyZBosVAt9Y6jTF3AquoO21xmbU2yxhze/3yxdbaHGPM28DXgJu6Uxs3t2XhIoFk6/4y0jOyWLv9ECN7deap68YycWBXf5clQcar89CttSuAFQ2eW9zg8ePA474rTSTwlTpqeWr1Nl7+dCedoiN48KrRzJ7Un3Bd5SltQFeKirQBt9vy+qYCHns7l0MVNVw/MZVfTRlO1/gof5cmQUyBLuJjX+05yryMLL7cc5RxqUm8dMtEzuiT6O+yJAQo0EV85FB5NY+9vYV/fL6HbvHR/OEHY7h6bB810ZJ2o0AXOU1Ol5u/rtvFk6u3Ulnj4rYLBvKLS4eSoCZa0s4U6CKnYV3+IdIzssgtKuPCod2Zd+UohvRUEy3xDwW6SCvsK6niocwc3vp6H32SYln8o/FcMTpZTbTErxToIqeg2uli6cc7WPR+Hm5rmXvpUH528WBiItVES/xPgS7ipfdz9zN/eTY7D1VyxehkfjtjFP26qomWdBwKdJEW7DxYwfy3snk/t5jBPeL5y60TuXBoD3+XJdKIAl2kGZU1Tha9n8fSj3cQFRHGb6aP5MfnDVATLemwFOgiDVhreevrfTy8Iod9JQ6+N7YPv542gp5qoiUdnAJdxENuUSnpGVmsyz/M6N6dWTR7LOP7q4mWBAYFughQUlnLH9/dyl/W7SIhJoKHrj6D6yakqomWBBQFuoQ0t9vyj417eGzVFo5W1nDDpP7815RhJMWpiZYEHgW6hKwvdh8hPSOLrwpKmDCgC+mzJjK6t5poSeBSoEvIOVBWzWNv5/L/Pi+gZ0I0T117Nled3VtXeUrAU6BLyKh1uXnl0108tXorDqeLn140iLsuGUqnaP0aSHDQT7KEhLV5B0lfnsXW/eVcNKwH9185isE9Ovm7LBGfUqBLUCs8WsVDmdms+KaIfl1j+dNNaVw2sqd2r0hQUqBLUHLUuvjTmnye/TAPgF9ePow53xmkJloS1BToElSstbybU8yDb2Wz+3Al089M4b7pI+nbRU20JPgp0CVo5B8o54Hl2Xy09QBDenbib7dN4vwh3f1dlki7UaBLwCuvdvLM+9tY9u8dxESE89sZdU20IsPVREtCiwJdApa1loyv9vLwihz2l1Zzzfi+3Dt1BD0Sov1dmohfKNAlIGXvrWuitX7nYc7qm8jzPxrPuNQu/i5LxK8U6BJQjlbW8Id3tvK3z3aRFBfFgu+dyQ/T+hGmJloiCnQJDC635X837OHxVbmUVNVy4zn9+eXlw0mMi/R3aSIdhleBboyZCiwEwoGl1toFzaw3AVgHXGutfd1nVUpI+3zXEeZlbGZzYSkTB3blgVmjGdmrs7/LEulwWgx0Y0w48CxwOVAAbDDGZFhrs5tY71FgVVsUKqGnuMzBgpW5vLGpkJTOMTx9/ViuPKuXrvIUaYY3n9AnAnnW2nwAY8xrwFVAdoP17gL+D5jg0wol5NS63Lz0yU4WvreNGqebn188mDsmDyFeTbRETsqb35A+wB6PxwXAJM8VjDF9gKuBSzhJoBtj5gBzAFJTU0+1VgkBH287QHpGFtsPVDB5eA/uv3I0A7vH+7sskYDgTaA3Nb+1DR4/BdxrrXWdbDpsrV0CLAFIS0tr+BoSwvYcruShzBzeziqif7c4/vzjNC4dmezvskQCijeBXgD083jcF9jbYJ004LX6MO8OTDfGOK21//JJlRK0HLUuFn+0nec/3E6YMdxzxXBuvWCgmmiJtII3gb4BGGqMGQgUAtcBsz1XsNYOPPa1MeYl4C2FuZyMtZZVWfv5fWY2BUeqmHFWL34zfSS9k2L9XZpIwGox0K21TmPMndSdvRIOLLPWZhljbq9fvriNa5Qgk1dczgPLs/h420GGJyfw6n9M4rzBaqIlcrq8Om3AWrsCWNHguSaD3Fp78+mXJcGozFHLM+/nsezfO4iNCmfelaO48Zz+RKiJlohP6DwwaXPWWv75RSGPrMzlYHk1Pxzfj3umDqd7JzXREvElBbq0qc2FJczLyOLzXUcY0zeRP92Uxtn9kvxdlkhQUqBLmzhSUcMT72zh1fW76RoXxWPfP4trxvdVEy2RNqRAF59yuS2vrt/NH97ZQpnDyc3nDeA/LxtGYqyaaIm0NQW6+MyGnYeZ92YW2ftKOXdQN9JnjWZ4SoK/yxIJGQp0OW37Sx08siKHf325l16JMSyaPZYZZ6qJlkh7U6BLq9U43bz4yQ6efm8btS7LnZOH8PPJg4mL0o+ViD/oN09a5aOtB3hgeRb5Byq4bGRPfjdzFP27qYmWiD8p0OWU7D5UyYOZ2azO3s/A7vG8eMsEJg/v6e+yRAQFunipqsbF8x/msXhNPhFhhv+eWtdEKzpCTbREOgoFupyUtZaVm4t4KDOHwqNVzBrTm/umjyQlMcbfpYlIAwp0ada2/WWkL8/ik7xDjEhJ4H/nnMOkQd38XZaINEOBLo2UOmpZ+O42Xl67k7iocOZfNZrZE1PVREukg1Ogy3Fut+WNLwpZsDKXQxXVXDehH7+aMpxuaqIlEhAU6ALANwUl3J+xmS92H2VsahLLbk7jrL5qoiUSSBToIe5QeTVPvLOF1zbsoVt8NE/8YAzfG9tHTbREApACPUQ5XW7+9lldE63KGhe3nj+QX1w2lM4xaqIlEqgU6CHos/xDzMvIIreojPOHdCP9ytEMTVYTLZFAp0APIUUlDh5ekUPGV3vpkxTL8zeMY+oZKWqiJRIkFOghoNrp4s//3sGi9/Nwui2/uHQoP7toMLFRuspTJJgo0IPcB7nFzH8rmx0HK5gyKpnfzRxFv65x/i5LRNqAAj1I7TpUwfzl2byXW8yg7vG8/JOJXDSsh7/LEpE2pEAPMpU1Tp79II8/rdlBZLjhf6aN4JbzBxIVoas8RYKdAj1IWGvJ/GYfD2XmsK/EwdVj+/DraSNI7qwmWiKhQoEeBLYUlZGekcWn+YcY1aszz1w/lrQBXf1dloi0MwV6ACupquWPq7fyl3W7SIiJ4PffPYPrJ6YSrqs8RUKSAj0Aud2W1z8v4NG3czlcWcPsian8aspwusRH+bs0EfEjrwLdGDMVWAiEA0uttQsaLL8BuLf+YTnwM2vtV74sVOp8ueco8zKy+GrPUdL6d+HlWRM5o0+iv8sSkQ6gxUA3xoQDzwKXAwXABmNMhrU222O1HcBF1tojxphpwBJgUlsUHKoOllfz2Nu5/GNjAT0SovnjtWP47tl9dJWniBznzSf0iUCetTYfwBjzGnAVcDzQrbVrPdZfB/T1ZZGhzOly88qnu/jju1upqnEx5zuDuOuSISSoiZaINOBNoPcB9ng8LuDkn75vBVY2tcAYMweYA5CamupliaFr7faDpGdksXV/ORcO7c68K0czpGcnf5clIh2UN4He1JzeNrmiMZOpC/QLmlpurV1C3e4Y0tLSmnwNgb1Hq3hoRQ6ZX++jb5dYXrhxPFNGJWv3ioiclDeBXgD083jcF9jbcCVjzFnAUmCatfaQb8oLLY5aF0s/zufZD7bjtpa7LxvGTy8aREykmmiJSMu8CfQNwFBjzECgELgOmO25gjEmFXgDuNFau9XnVYaA93L288DybHYfrmTq6BR+M2OkmmiJyClpMdCttU5jzJ3AKupOW1xmrc0yxtxev3wxcD/QDXiufreA01qb1nZlB48dByuYvzyLD7YcYEjPTvz11klcMLS7v8sSkQDk1Xno1toVwIoGzy32+Po24DbflhbcKqqdLPogjz9/vIOoiDB+O2MkPz5vAJHhaqIlIq2jK0XbmbWWjK/28siKXIpKHXx/XF/unTacnglqoiUip0eB3o5y9pUyLyOL9TsOc0afzjx7wzjG9+/i77JEJEgo0NtBSWUtT67ewl/W7SIxNpKHrz6Tayf0UxMtEfEpBXobcrkt/9i4h8dXbeFoZQ0/Oqc/v7x8GElxaqIlIr6nQG8jm3YfYd6bWXxTWMLEAV1JnzWaUb07+7ssEQliCnQfKy5z8OjKLfzfpgKSO0ez8LqzmTWmt67yFJE2p0D3kVqXm5fX7mThu9twOF3cftFg7rpkCPHR+i8WkfahtPGBT/IOMi8ji7zici4e3oP7Z45iUA810RKR9qVAPw0FRyp5KDOHlZuLSO0ax9Kb0rh0ZE/tXhERv1Cgt4Kj1sULH+Xz/Ed5APzX5cP4j++oiZaI+JcC/RRYa1mdvZ8HM7PZc7iKGWf24r4ZI+mTFOvv0kREFOje2n6gnAeWZ7Nm6wGGJXfi1dsmcd4QNdESkY5Dgd6C8monz7y3jWWf7CAmIpz7Z47ixnP7q4mWiHQ4CvRmWGv515eFPLIil+Kyan4wvi//PXUEPRKi/V2aiEiTFOhNyNpbwrw3s9i46whj+ibywo3jGZuqJloi0rEp0D0cqajhD6u38Opnu0mKi+LR75/JD8b3I0xNtEQkACjQqWui9ff1u3ninS2UOZzcdO4A7r5sGIlxkf4uTUTEayEf6Bt3HmZeRhZZe0uZNLArD1w1mhEpaqIlIoEnZAO9uNTBIytz+ecXhfRKjOGZ68cy86xeuspTRAJWyAV6jdPNS2t3sPDdbdS6LHdMHswdk4cQFxVy/xUiEmRCKsXWbD1A+vIs8g9UcMmIntw/cxQDusf7uywREZ8IiUDfc7iSB9/K5p3s/QzoFseym9O4ZESyv8sSEfGpoA70qhoXz3+0nRc+2k6YMdxzxXBuu3Ag0RFqoiUiwScoA91ay6qsIh58K4fCo1VcOaY3900fQa9ENdESkeAVdIGeV1xGekY2/+51FBUAAAZSSURBVM47yIiUBP7+H+dw7uBu/i5LRKTNBU2glzlqWfjuNl5au5O4qHDSrxzFj87pT4SaaIlIiAj4QHe7Lf/8opBHVuZyqKKaa9P6cc8Vw+nWSU20RCS0eBXoxpipwEIgHFhqrV3QYLmpXz4dqARuttZu8nGtZOZnsmDzk3QaUcxPP0zm2sE/5Z31vdm0+yhn90vizz9OY0y/JF+/rYiIT2TmZxI1YAF/Ky7h/ddTmDtuLjMGzfDZ67cY6MaYcOBZ4HKgANhgjMmw1mZ7rDYNGFr/ZxLwfP3fPpOZn0n62nQcLgfGwEHHfhZ98wgRjh/y2DU3cM24vmqiJSId1rEMM5EOAPZV7CN9bTqAz0Ldm0/oE4E8a20+gDHmNeAqwDPQrwJesdZaYJ0xJskY08tau88nVQILNy3E4XKc8JwJq8WduII/rTmXP63J99VbiYj4XHHSo7jDT8wwh8vBwk0L2zXQ+wB7PB4X0PjTd1Pr9AFOCHRjzBxgDkBqauopFVpUUdTk8zbiCEOTO53Sa4mItLci55Gmn28m21rDm0Bvaj+GbcU6WGuXAEsA0tLSGi0/mZT4FPZVNP7A3yu+F89dM/5UXkpEpN1Neb1XkxmWEp/is/fw5py+AqCfx+O+wN5WrHNa5o6bS0x4zAnPxYTHMHfcXF++jYhIm2iPDPPmE/oGYKgxZiBQCFwHzG6wTgZwZ/3+9UlAiS/3n8O3Bw0WblpIUUURKfG+P0IsItJW2iPDTN1xzBZWMmY68BR1py0us9Y+ZIy5HcBau7j+tMVFwFTqTlu8xVq78WSvmZaWZjduPOkqIiLSgDHmc2ttWlPLvDoP3Vq7AljR4LnFHl9b4I7TKVJERE6ProsXEQkSCnQRkSChQBcRCRIKdBGRIOHVWS5t8sbGHAB2tfKfdwcO+rCcQKBtDg3a5tBwOtvc31rbo6kFfgv002GM2djcaTvBStscGrTNoaGttlm7XEREgoQCXUQkSARqoC/xdwF+oG0ODdrm0NAm2xyQ+9BFRKSxQP2ELiIiDSjQRUSCRIcOdGPMVGPMFmNMnjHm100sN8aYp+uXf22MGeePOn3Ji22+oX5bvzbGrDXGjPFHnb7U0jZ7rDfBGOMyxlzTnvW1BW+22RhzsTHmS2NMljHmo/au0de8+NlONMYsN8Z8Vb/Nt/ijTl8xxiwzxhQbYzY3s9z3+WWt7ZB/qGvVux0YBEQBXwGjGqwzHVhJ3R2TzgE+83fd7bDN5wFd6r+eFgrb7LHe+9R1/bzG33W3w/c5ibr79qbWP+7p77rbYZvvAx6t/7oHcBiI8nftp7HN3wHGAZubWe7z/OrIn9CP35zaWlsDHLs5tafjN6e21q4Dkowxvdq7UB9qcZuttWuttcduTriOurtDBTJvvs8AdwH/BxS3Z3FtxJttng28Ya3dDWCtDfTt9mabLZBQf3+FTtQFurN9y/Qda+0a6rahOT7Pr44c6M3dePpU1wkkp7o9t1I3wgeyFrfZGNMHuBpYTHDw5vs8DOhijPnQGPO5MeamdquubXizzYuAkdTdvvIbYK611t0+5fmFz/PLqxtc+InPbk4dQLzeHmPMZOoC/YI2rajtebPNTwH3WmtddR/eAp432xwBjAcuBWKBT40x66y1W9u6uDbizTZfAXwJXAIMBlYbYz621pa2dXF+4vP86siB3iFuTt3OvNoeY8xZwFJgmrX2UDvV1la82eY04LX6MO8OTDfGOK21/2qfEn3O25/tg9baCqDCGLMGGAMEaqB7s823AAts3Q7mPGPMDmAEsL59Smx3Ps+vjrzL5fjNqY0xUdTdnDqjwToZwE31R4vPoQ1uTt3OWtxmY0wq8AZwYwB/WvPU4jZbawdaawdYawcArwM/D+AwB+9+tt8ELjTGRBhj4qi7+XpOO9fpS95s827qZiQYY5KB4UB+u1bZvnyeXx32E7q11mmMuRNYxbc3p87yvDk1dWc8TAfyqL85tb/q9QUvt/l+oBvwXP0nVqcN4E51Xm5zUPFmm621OcaYt4GvATew1Frb5OlvgcDL7/ODwEvGmG+o2x1xr7U2YNvqGmP+DlwMdDfGFADzgEhou/zSpf8iIkGiI+9yERGRU6BAFxEJEgp0EZEgoUAXEQkSCnQRkSChQBcRCRIKdBGRIPH/AXclHM/6R8xkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.triplot(points[:,0], points[:,1], tri.simplices)\n",
    "plt.plot(points[:,0], points[:,1], 'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59cc5e85",
   "metadata": {},
   "source": [
    "### in 3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "a08fb110",
   "metadata": {},
   "outputs": [],
   "source": [
    "# u,v parameterisation variables\n",
    "u = np.array([0,0,0.5,1,1]) \n",
    "v = np.array([0,1,0.5,0,1]) \n",
    "\n",
    "x = u\n",
    "y = v\n",
    "z = np.array([0,0,1,0,0])\n",
    "\n",
    "tri = Delaunay(np.array([u,v]).T)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b4baea53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3xb9bn/P1qWJXnEI45X7MQjXpmOnTi0UBpGb3PbtJcu6KW5KQ2UEZqUwiWFltkCofNXoIX2cqG0DdBfW5oCGcxQaIhXFlkekmzZ8pRk7XnG74/8zkGSNY6kI8myz/v1yotgSed8Tyx99Jzn+3yeR0TTNAQEBAQEUoM43QsQEBAQWEgIoisgICCQQgTRFRAQEEghgugKCAgIpBBBdAUEBARSiDTK40Jpg4CAgEDsiMI9IES6AgICAilEEF0BAQGBFCKIroCAgEAKEURXQEBAIIUIoisgICCQQgTRFRAQEEghgugKCAgIpBBBdAUEBARSiCC6AgICAilEEF0BAQGBFCKIroCAgEAKEURXQEBAIIUIoisgICCQQqJ1GRMQCAtN06AoCh6PBwRBQCqVQiwWQyKRQCwWQywWQyQK22xJQGBBIooymFJo7SgwC5qmQZIkCIII+DvzmL/QMiLM/BHEWGCBEPYNLoiuAGeCxVYkEsHtdkOr1YIgCOTk5ECpVEKlUiErK4t9DfNneHgYJSUlUCqVghgLzHfCvpGF9IJAVGiaBkEQGB0dRV5eHlQqFZxOJzQaDZxOJ6qqqiCRSOByuWA0GqHT6eD1eiGRSKBSqaBSqaBUKuFyuUDTNMTii1sJBEHA5/MFnEsQY4H5jhDpCoSFEVsmddDX14ecnBwYDAYQBIHly5ejqKiIjX4ZMWUgCAJOpxMOhwMOhwOTk5MQi8WQyWQBYqxSqZCdnR1wXuYPI7YikQgSiYTNGzPiLIixwBxFiHQFuENRVECeFgDMZjOmpqZgsVjQ2NiIRYsWsY+FEz6pVIq8vDzk5eUBuCjC5eXlbKTscDhgsVgwNjYGt9sNsVjMijAjyAqFgj1ecGqDiZrDRcaCIAvMRQTRFWChKAoEQYAkSfZnBoMBWq0W2dnZKCoqwpIlSwIENxYYoZRIJMjNzUVubm7A4yRJwuVyweFwwGq1Ynx8HG63GwBCijETWQd/QTDnstlsKCoqEsRYYE4hiO4Ch7mN9/l8oCiK/dnk5CSGh4eRl5eHVatWQalUYnBwEFHSURFhRDccEokEOTk5yMnJCfg5RVGsGNvtdkxOTsLlcgEAFApFgBgrlUpWVPv6+rB+/fpZa2CiY/9UhSDGAqlCEN0FClNjSxBEgNiOjY1Bp9OhqKgI69atC8i1ikQi9rnxEE10wyEWi1lh9YeiKLjdbjZnPD09DafTCZqmoVAo4PV6YTAYoFKpoFAoIJFI2OsmSRJerzfgeP5pCiY6FsRYgG8E0V1gBIutSCQCSZIYHR2FXq9HaWkp2tvb2ZIvf+IVzWTB5ICVSiUWL17M/pymabhcLpw6dQputxtGoxFOpxMURSE7O5sVcCY6Dhbj4FrjUDljoaJCIF4E0V0ghKqx9fl8GBkZwcTEBCorK9HR0QGpNPxbQiwWJzW9wBcikQgKhQJSqRTV1dXsz2mahtvtZjfxZmZm4HQ6QZIk5HL5LDGWSqURxVgobxOIB0F05zmM2JpMJthsNlRUVMDr9WJoaAhGoxFVVVW45JJLZpV7hSJd6QW+YMRYoVCgqKiI/TlN0/B6vbDb7XA6nRgbG4PD4QBJksjKypolxjKZLCAX7vV6IRKJoFarUVtbK4ixQEQE0Z2nMDW2TIRGEARMJhMsFgusViuWLVuGFStWcBJbBj5Ecy6lJxhEIhHkcjnkcnlIMXY4HHA6nZiYmIDD4QBBEAG1xswfi8UiGD8EoiKI7jwj2NAgEongcDig0Whgt9vR0tKC5ubmuD7oYrE4ZKTL9VipFpdEz+cvxoWFhQGP+Yvx1NQUu5l34sSJWWIsk8nY14UTY7FYLBg/FgiC6M4Tgg0NIpEIFosFGo0GFEWhtLQUZrMZJSUlcZ8jk9ILyT5PVlYWsrKyUFBQwJ6vp6cHq1atYgXYYDBgeHh4liWa+eO/WSkYPxYOguhmOKEMDSaTCVqtFlKpFLW1tcjPz4fNZsPMzExC50pUNNOd000FMpkMixYtmmUgIQiCFWOTyYSRkRF4PB5IJJIA44dKpYJcLmdfF874EUqQBTHODATRzUBCGRoAYHp6GlqtFiqVCk1NTQEmg3CpgVhItHphPhNc2RCMVCpFfn4+8vPzA35OkiQrxmazGXq9PqQlOrg/BUVR0Gq1UCgUbLmcYPzIDATRzSDCGRomJiYwPDyMRYsWYc2aNQH9Chj4EN1MSi+kmmiiGw6JRBLQn4KBJMmA/hTj4+NwuVwQiUSsGDscDshkMlZYBeNHZiCIbgYQytBAURT0ej1GR0exePFirF+/PuC2NBi+RFdIL4QmXtENR7j+FBRFsWJsNBoxOTmJ8fFxAIGWaMaFx9ydBNcaM/8VjB+pRxDdOQxTYzs9PQ0AKCgoAEmSGBkZwdjYGMrKyrBhw4aA3fFwCKKbXPgW3XCIxWK2P4XVakVxcTEKCgoC+lP4W6IBxOXCE8rbkocgunOQYPeY3W4HSZIwGAyYnp5GZWUlNm3aBIlEwvmYfOV0wx1joX8YUyW6/lAUxdYFc+1PYTAY4HK5olqig40fDG63G3K5nO1lIYhx7AiiO4cINjSIRCJ4PB5MTk7CbrdjxYoVqKuri8nQwLAQI91UCkG6RTcckfpT+IuxyWQK6E8R3EqTsUSPjY2hsLBw1rUKkTF3BNGdA4QyNLhcLmg0GjgcDixatAiFhYWoqKiI+xx8vPnDVS84nU6YTCbk5uYG9LkNtYb5Uqcb6nxzUXTD4W+JLi4uZn9O0zQ8Hg8rxnq9PqA/BVNZIZVKoVKpWDEGohs/BDG+iCC6aSSUocFms0Gj0cDr9aKmpgZFRUWYnp6GxWJJ82pnVy8wTjfmi4Hp5gV8vKmTk5PDbuoAc9MGzAeZJrrhEIlEyM7OZpvWMzCW6HPnzoGmaYyPj7P9KWQyWcBQUsaFx/yuIxk//EvbFkpFhSC6aSDY0CASiTAzMwONRgORSISamhrW6QRc3MlONDXAB8wHxm63Q6PRwOVysV8MPp+P/bD4b+rYbDZMTEzA5XKxPQtIkgyoPZ0PH7L5IrrhYCzRUqkUFRUVAV+iPp8vYA5ecH8KfzHOysqaJcbB5wkXGc+H9wkgiG7KCGdoMBqN0Gg0yM7ORkNDw6wSIeDiLZq/4yxduFwumEwmuFwu1NbWsrm94OjVf1PH33as1+vhcrkgl8tnzUbz39DJyclBVlZWRn3I5rvoMpAkGbCBKxKJZlmiGfzFeHp6GkNDQ/D5fGxqIliMGSiKgsfjCTiOxWLBkiVL5oXxQxDdJBPK0AAAk5OTGBoaQm5uLlauXDlr19mfdEe6NpsNg4OD8Hg8UCgU2LBhQ1zHkUgkyMrKQmlpacDP/V1ZMzMzGB0dZS2yocR4LrJQRTcS4SzRPp8voNZYp9PN6k/BCLJcLgdBEDAajVi8eHFU4wcT1ASbTeYSgugmCWZzzGq1snO7mFyYTqdDYWHhrHE44UhXpGuxWKBWq0FRFGprayGXy9HX18f7ecK5svz7FRiNRgwPDwdESky+OLiTVzpIl+hmotDLZLKQlmiCIAIazDNfvsxd4ujoKCvGTFrKv9YYAA4fPoxTp07hkUceSWiNyUQQXZ7xr7ElCAJnzpzB+vXrodfrodfrUVJSgra2tpgiNj7KvWLBbDZDrVYDAGpra9lIxe12p7RkLFy/gkg5RP+Nu1Ru2jGbQ6kmHedMltBLpdKQX77T09MwGAzIysqalZZiyuFUKhWcTidmZmZmvV/mGoLo8kSocTjMrK6uri6Ul5dj48aNEcfhhEMikfAW6UaKyGZmZqBWqyEWi1FXVzfrzTtXei+Eu21lpj84HA52d72rq2vWKB6VShWTsYQL87UqYy5AURSUSmXItJT/hu3evXtx7NgxUBSF06dPo7m5GTt37uQkwjfccANee+01lJSU4MyZM7Mep2kau3btwoEDB6BUKvH888+jtbU1rusRRDdBQhkafD4fhoaGYDAYIBKJ0NHRkdCHnK+cLhMxB6/FZDJBrVZDKpVixYoVYfNhc90ckZWVhcLCQhQWFoIgCLjdbqxbty5AjPV6PRwOxyxHFlPyFG/kmI5b/XSQji8XgiBCBisSiYS1RAPAM888g71796K5uRlr1qzBuXPnON9Rbt++HTt37sS2bdtCPn7w4EEMDAxgYGAAnZ2duOWWW9DZ2RnX9QiiGyehDA1utxtarRZWqxXV1dWor6/HsWPHEo6q+Mrp+osuTdOs2GZlZaGxsTFk5UTw6zMtoos0isffkcXUGDPj24NrjLnOkEslmfa7iBeCICI2c/LHYrGguLgYLS0taGlp4XyOyy67DENDQ2Ef379/P7Zt28YGUWazGePj4ygrK+N8DgZBdGMklKGBMQm4XC4sX7487nE44eBL7BjRNRgMUKvVUCgUaG5uDui7G4m5kl7gg0iOLOaW1W63Y2pqCi6XCwBm9bdVKBTs7zkdG2mpJl3RPEmSnNNyVqt1VtqJD/R6PZYuXcr+f2VlJfR6vSC6ySSUoYEZh0OSJGpqakJ60ucKjKOot7eXU5laKOZ6eoEPmH61wb0KGMOH3W4PMHwwmzkikQgkSbINYVLxPkj1ey2WcjE+CZdeCIXVak3KRlqo9228//6C6EYgnKGBcY/5j8MJBxMdpmOXGbh4DVNTU9BoNPD5fFi5cmXAbXYszNdpwFwI18WLaTY+MTEBp9OJ/v7+WYYPJk3Bp+EjHf+OofYDUkEsomuxWGaZNPigsrISIyMj7P+Pjo6ivLw8rmMJohuCcIYGxlWjUChmjcMJB1N5kGrRpWkak5OT0Gq1yMvLw9q1a9Hf3885NxYKPqbrzjeYZuNutxtSqRTLly8HEGj4MJlMAQYA//riYDcWV9KRzkhnpMv1vBaLJSnpha1bt+LJJ5/Etddei87OTuTn58eVWgAE0Q2AKfti6j+rqqpY8RoaGsKiRYuwatUqKJVKzseUSCRsDWkqYMb3aLVaLFq0KMCAkep632AyIb0QL8EiGM3wYbfbA6yxTJ+CcKPbg0lHfjVdohtLTpdpwBMr1113HY4cOQKDwYDKyko8+OCDbMe0m2++GVu2bMGBAwdQV1cHpVKJ5557LuZzMAiii9k1tiRJwmg0QiQSYWRkBMXFxVHH4YQjVTW2FEVhYmICQ0NDKCgoQGtr6yy320IS3VRHglxTSJEMH0xZm7/hIysra1ZZG9NOcS5bgPmEa3ohkffWiy++GPFxkUiEp556Ku7j+7OgRTeUoYGiKIyNjcFoNCI/P5/zOJxwSKVSXkTXf66VP8x6h4eHo345pFt05zOJirxMJkNBQUFAPpLZ/AzV21Ymk8HlcmFiYiJg6kMySUeajIHLvy0junM9jbUgRTdUjS1BEBgeHsbU1BTKy8uRk5OD2trahM/FV6QbnBtmBlPqdDosXrwY7e3tUXODyRRdLm/0ZEe66d6w5PvD7l9jXFhYGHAupi+G1+sNmPoQPJwyEcNHMOnaSOOK2+2OKfWXLhaU6IYSW4/HA61Wi5mZGVRVVWHTpk0QiUTshNVE4VN0KYoCRVEYHR3FyMgISkpKOIktQ7oj3WSK7rHtP4Gh8zw+ffgxqKpKor+AZ1I9hkgmkyE7OxtVVVUBa3C73WyawmAwzGoqHzwpOBbSlV7gitlsntPdxRgWhOiGMjQ4nU5otVrYbDYsX74cjY2NSbkt4Ut0RSIRdDodpqamsGTJkrjSHvNVdD/cthf6/R8CAA6vvxWX/v1BLNqwgvfzRCPdOWR/w0eoGmMmTRHK8MFUVERqKp8O0Y1lwzBZlQt8M69FN5ShwX8czvLly9HS0pLUD0uiosuMXDeZTFAoFHE3zQGSJ7perxd6vR7Z2dnIyclJmTkAuBjhMoILAJSXwHv//gOs+emNwPr4SnriIdUbW7GkUiJNCmZaKVqtVoyPjwcYPvxL2+RyOUiSTHk/41hrdOd6hzFgHopuKEODSCSC2WyGRqMBgFnjcCIdK1HxiFd0CYLAyMgI9Ho9KioqUFJSgvLy8rgFF0jcxhsM09hnamoKJSUlsFgs0Ov1bANy5kPL/JfvSLfzhp9h9JV/zX6ApnHqe79F3hfbgd+v4+18kZir1RKREIvFAQ1jGBjDh91uD+hr6/P5oFQq4fP5AmqMk3ndsdToms1mIdJNJeEMDUajEVqtFllZWaivr+ec8wnXkStWJBJJwOiRaDAbeuPj46isrMSmTZsgkUhw7ty5hAWTr0jX5/NheHgYk5OTqK6uRkdHR0DqhrkOph6VKYHyer0gCAKDg4MBYhyPeHTd9AuM/PX9iM+x/r0b7/37D3DZaw8nXRAzUXTDwRg+ghsgDQwMIDs7GxKJJGDiA9NU3r+0ja+69FhqdIX0Qopgyr5IksTJkyexevVqAMDU1BSGhoaQk5ODlpaWmPsMSKXSmL5lIx3H4XBEfZ6/kPmLLQMfueFEu5URBAGPx4POzk5205H54DOtLRlC1aPa7Xao1WoUFBTA4XBgZGSE3XX3zy3m5OREzC123/x/oHv5PU5rnv7gDA6u/jau/OAXyMqP7T0QC6munEhHpQZN08jNzQ05fsd/FppWqw0YTOn/BRvrnVqsfRcE0U0ioWps3W43RkdHMTo6ioKCAqxdu5bTOJxQ8F3qFQ6v18uWqi1dujRAyPzhI0oVi8WsyyYWCIKATqfD+Pg429ounjQHM3a7qKgooP9D8PTg8fHxgP4F/mmKU7ufxvCL78Z0XqduCgead+DTb+9FfmNV9BfEQToi3bniSIvUVJ4R44mJCdjtdjYvHGyFDhfcxJrTjdeam0oyTnTDTWgYGRmB3W6H3W6PeRxOKJhIN1HCia7X68XQ0BCmp6dRXV0dVmyjHScWYhVukiSh0+kwNjaGiooKdHR0oKenJ+7zh8vphpsezFiy7XY7DAYDum97ArYDJ+I6N2F34a1LdmPTn76P8s+2x30N4ZhP6YVwxFq9EGpKsL/hw263h20qz9QYxyq6QqSbBEiShM/nY1vp6XQ6TExMoLy8HMXFxVi6dCkvO6zJEl2Px4OhoSEYjUZUV1ejrq6O04eHr0iXyzGYigm9Xo/y8vKAyReJbIbF+lr//gW9u34dt+Ay0CSFo9f+GCsf+AYav/ulhI4169gpFt1MtQFHMnz4N5VnDB8ejwcymYzdvMvJyQlbYyyIbpJgbpGZcThLly5lReH8+fNx3T6Hgm/RZaZKzMzMYNmyZaivr4/pQ5OKSNffeFFWVhayPC0dTWuO3/E0tM+/wdvxzjzwB1jODGPjs3fwdsyFEukm65zhmsprtVpIpVJkZ2dHNHxMTU0JopssLBYLTp48GVK4+BJKgL+crs/ng8ViwfHjxxMyYUgkEni93oTWEk50/S3FpaWlEWuBE4m44xHsE3c+A82zh+I6XyRG/vJPTH80iDV/vht5BYvi2uTxJx2im8h64z1nqs0RJEkiPz8fhYWFYQ0fU1NTuPfee3H+/Hl8/etfx8qVK3HZZZfhxhtv5HSOQ4cOYdeuXSBJEjt27MCePXsCHrdYLLj++uuh0+lAEATuvPNOfPOb34z7mjJOdPPy8lirbjB8im6ix3K5XNBoNLBYLJBIJGHXzBU+5qQFC6Z/s5ySkhJOLrdUphdO/vfvoP7dwbjOxQV33xh6N/8AtS/cDq/i4pcsk1dkNnq49i5IdfSfruqFdETXob5cgvcB9u/fj82bN+Ott97C8PAwZmZmOB//tttuw5tvvonKykq0t7dj69ataG5uZp/z1FNPobm5Ga+++iqmp6fR0NCA//zP/4w7jZlxoisWi8OKF5P74QOpVBpTfS2D0+mERqOB3W5HTU0NmpqacOzYMV5MFnzldCmKwvj4OIaGhjg3y/E/RipE9+T3n8XgM6/HdZ5YIIw2DHzpp/jU6w+jsH0Fm1dkNu+YW1l/h1Y41918LxlLB7GUbfp8PuTl5bFlo1zo6upCXV0dampqAADXXnst9u/fHyC6jJOVpmnY7XYUFhYmdJeRcaIbiXiFMtyxuNTXMjDDKR0OB2pra3m3F/MR6TJDNI8dO4aioqKYxNb/GKHEn89c76l7n8Pgr1/l5VhcoNxevHvVHrT9eheWXXf5rLwiY5e12+1hXXdut5u3vsmc1pymSDfVJLuXbqiBk8Gj1Xfu3ImtW7eivLwcNpsNL7/8ckL/9hknupGEjM9Il5n4EA273c5OAq6trUVRUVFScnuJ5FKZaRJqtRoURWHjxo1xj+1Jdnrh9A+ex8CT++M6fkJQNHpu/iWs54ex+qH/CngonF2WIAi2o5fT6YRarcbg4CDkcnlAiiJe113E5S6gSDeWqDLWzx6XgZOHDx/G2rVr8c4770CtVuOqq67CpZdeGndHs4wTXSD8h5fvnG6kyMVms7H9TGtra5M+CTiejT3/OWnMqCG1Wp3QnLRkphc+uv8F9D/x93iXxgv9/+cVWM4O45N/+WHU36dUKmVNARaLBVVVVVCpVPB6vWzNuMlkYu+YmN12Lq67aKRadNM1Yp7rdTJTmGOFy8DJ5557Dnv27IFIJEJdXR2WL1+OCxcuYMOGDTGfD8hQ0Q1HKjbSrFYr1Go1CIJgxTYVxBLp0jSN6elpqNVq5Ofns3PSPB5PwnnhRJvmhBPdMw/9EX2//Fvcx+WTybeO4/D623Dl+z+DVKXg9BpGlPzrUEO57pgR7qFcd/5Tg6OxEGzHDFzEPt4OY+3t7RgYGIBWq0VFRQVeeukl7Nu3L+A5VVVVePvtt3HppZdicnISfX19bA44HuaV6PK9keYvukynfoqiUFtbG9OYZz7GsHOJdGmahsFggFqtRm5uLtauXQuF4mPR4MNgkWh6IRRnfvQnXPjZXxJZFu/Y1WN4vWkHrnjvZ8hZXhr1+dEiwUgj3MMNqvQX4mCrbKptwHO9gXm8NbpSqRRPPvkkPvOZz4AkSdxwww1oaWnB008/DeDiUMof/vCH2L59O1atWgWaprF3796AnH/M54z7lWkkFekFRuTMZjPUajUAoLa2Nq5fLB9j2CNVL9A0DaPRCLVaDZVKhdWrV4ccWzIXRDf4tWcffQkXfvJ/E1pTsvBZHDjcvhOf+L8/QOmn10Z8bry33+GmBvunKPytskyKwu12w+PxQKlUpkR80yG6sbzPrFZr3L10t2zZgi1btgT87Oabb2b/Xl5ejjfe4M+ck5GiG45E8o3BWK1WWK1WaLVa1NXVJdQcmRHdRNrdhapeoGkaJpMJg4ODUCgUUcfDp9JKHIpgcTj3+Ms4/9hLCa0n2dA+Ah/8x4NY89gNqL/58+Gfx3POMysrC4WFhbOssowhYHJyEiMjIxgcHAxoOu6fouBzPekYShnL3aHZbM6IBuZAhopuMr/ZTSYT1Go1O4Nq3brEm2Anw8LLiK1cLsfKlSs5ta7k49+Nr9Kw8z/9C879OPLY6zkDTePU3c/CcnYYbU/sDPOU5G80iUQiKJVKKJVKjI6Oorm5GTKZLKDpeHCf2+BG8vHWl6bDjTYfm90AGSq6fMNEjMzOfmNjI3Jzc3H06FFejs+X6NI0jZmZGQwODkImk6G5uXlWGVOyCSe6sQjO+Z/9BWcf/iOfy0oJQy+8BeuFEVx+6BGIgwQonb0XwjUd9/l8bEnb+Pg4HA7HLNddpAYy/qQjvRCL6GbK1AggQ0U32oYF1zcIs/Gk0WigUCiSJmJ8iK7ZbIbD4cDQ0BD7pZAOgiNu/0oJiqLYDzLzJzilYnn5Qwz/7p1UL5s3TF19OLjyJlz5wS8hL/r4dzAXG97IZDIUFBTMaq3o77qbnp4OGFIZznU310XXarWirq4uySvih4wU3UhwmfjACIVGo4FKpQp7e85H1QHA3WgRCovFgsHBQYhEIt7SHYngH+kajUZ29M6qVasgFovZ21z/CQJyuRw5OTmYefkozBksuAyuMSMOrNyBzW/tRX7LMgDpqWON53zhunlFc90x1RJMZUUqiEXohfRCGmHKxkIVStM0jampKWg0GuTm5obd5WdgBJyPhuixRrpWqxWDg4OgaZrdyOMr3ZEIjJW4u7sbMpmM/cIiCAIEQczaiWeaVp9/4hUM7f1rGlfOL6TTg7cuvQMbn7sTlV+4JG3mAb6I5rrT6/Vwu9346KOP2C/SZLvu5uOoHiBDRTfSmztU2Zi/Mys/P39W/Wq0YyUqurGkF2w2GwYHB0GSJOrq6ubUG8lms0Gv10MkEmHVqlWcbJAikQi6597A4MNzu0ohHmiSwrFtj6Npz7Wgr6rLaNENB+O6s9vtKCgoQHl5OWiahsfjYVMUwa47/xRFIq67+Th+HchQ0Y2EVCplDRI0TbPdtAoKClhnFldSNScNuNjDYXBwED6fD3V1dWHNF+mIqJxOJwYHB+F2u1FcXAyVSsXZdz7429dx6u5nk7zC9HL+sZeg+Fcj8I/WdC8lafiXPDKpruzsbE6uO4lEMmtAJZdAhhluyQWLxRKTYSmdZKTocml6o9frMTw8jMLCQqxfvz4uXzaf0yPCdT9zOBwYHByEx+NBXV1dRFsxX2PhuQq32+2GWq2GzWZDXV0dioqKoNfrOdfpqv/nIE7e9buE1popuN6/gLc+sRtXHPkpJPLEx0XNNbjkV8O57giCCJgWzNV1R5IkpztSQBDdtEFRFOx2O3Q6HcrKyuIWW4ZkDqdkulI5nU5WbKMJIeNKS0R0uQi31+uFVquF0WhEbW0tmpub2bVx7b2gfu4wTnzvmbjXmYlYz+nwetO3cOU/fw5l5eLoL8ggEqlekEqlyM/Pn3X7H8l1l5OTA6vViuzsbE5Bgsfj4SzQ6WZeiK7/uBmFQoGKigrU19cnfNxkiC7T5JzpuxtLK0imHC5RZ1s40SUIAkNDQ5icnMSyZcuwYsWKWWvjUqereeFNnNj9m7jXmMl4jTYcar0Vl+1/EMWbmqO/IA7S0dc2GSVjkVx3TAXMxMQERkZGIrru0vHvkQgZKbrMB1z0g6AAACAASURBVNx/kCIzbsZsNnMe1RENPnO6Ho8HZ8+ehc1mQ21tLYqLi2POzSbLxstMVdbr9aiqqoo4Dj6a1Vr7x7dw/PanElpjpkN5fDjy2XvR+qtbUbPtKt6Pn468fqpswP6uu8nJSdTW1kKpVIZ13Q0MDOC9994DQRDo6upCS0tLTLX20eajAcCRI0ewe/du+Hw+FBcX47333kvoGjNSdAFAp9NhZGQEpaWlAbO9ZDIZr+0dE+1a5na7odPpYDAYsGrVqoBb9Vjh207sf4dQVlaGTZs2RY1mIqUXhl58F8dvezKh9c0baBrHb38K1jNDWPs4twGJXElHm8V0DaVkqhfCue7q6uqQl5eHEydO4H//939x9uxZ/Nd//RenoZRc5qOZzWbceuutOHToEKqqqjA1NZXwdWWs6Mrl8pBTa/kQSv9jMW6dWPF4PNBoNDCbzaioqABJkigpKUloPXxFuiRJYmxsjJ2RxmUgJUO427mhF9/B8VufSGht85HBZ16H5bwOl/3jId6i03SNX5+LjrTCwkJ0dHSgoqICzzwT2x4Cl/lo+/btwzXXXIOqqioASPgzDAAZO++jrKws5C8k3ROBvV4vLly4gN7eXixatAgdHR1YsmRJwmIJJB7pMkaFkydPwmKxYP369aivr48pRxxKdIdeehedN/0i7nXNd6b/+REOrrkZPpuTl+MtFNHlOn04XjdaqPloer0+4Dn9/f2YmZnB5ZdfjvXr1+OFF16I+TzBZGykGw6+56RxFTn/Hf9ly5ahoaGBjWz4yg0nEukyll2v14uGhoa4v7GD1zD05/dw9IafAJm1l5FynMOTeL3pW9j8zk+Qt6IyoWOlK70wV2eyxWuM4NK4iSAI9Pb24u2334bL5cKmTZvQ0dGBFStWxL3ejBXdcLdqfNyCM3CJdH0+H4aGhjA1NYVly5ahvr5+1pszlSaLYMxmMwYGBljL7vDwcELVD/6R7vBf38fR7Y8LgssRwubCGx3fQeNT38bSLR1xt1pMlwDOVcddvB3GuMxHq6ysZA1BKpUKl112GU6dOrUwRTccfL4xIomuz+fD8PAwJicnUV1dHXHHn681RZoeEQxjJ6YoCg0NDayDLNESG+b1ulc+wL++8RiQYeU6aYekcOHm38B8+xAU/9HGtlr078ymUCgivmdSPaonHcTyHo3XGMFlPtoXvvAF7Ny5EwRBwOv1orOzE9/97ndjPpc/8050+SSU6BIEgeHhYUxMTGDp0qURxZZvQk2PCMbfsltfXz/rzZjonYBYLMbM26dx8u4/CIKbABNPHETVlBPtz+yG2+1mTQKTk5NwuVysddZfjJmomKbpOT2vjA9ibXYTzzh0LvPRmpqa8G//9m9YvXo1xGIxduzYgZUrV8Z8roDzJvTqNBLpm57PloyMyBEEAZ1Oh/HxcVRWVqZUbBkiCWYoy26of6NERXfycC+G735BSCnwgO7l92DtG8Xmtx6DYvFiLF78sYuNsc4yQqxWq9momGkV6nQ6o0bFmUq09qz+WCwWtrogVqLNRwOAu+66C3fddVdcxw9FxopuJPjqDsbcSmu1Wuj1elRWVqKjoyNtUUaoHg7+G3g1NTVR64ATEd3R1zvRs/1nguDyiPmkGq+33IirPvgFsks+visJZZ1lGpCPj4/DbDZDo9HA6XTOGuMeqnl8IqQjneFfoxsNs9mM1atXJ3lF/JGxosulvWMiokuSJEZGRtiWdVyMA9FI1EnkL5hcLLvRjhEL+sPd+OdXHxJSCknAM2nGwVXfxqcO/hiFreHt60wD8tzcXIhEIixfvhzA7DHuwc3jYxnLE4p0DKWcr710gQwW3UgkUjZGURRGRkYwOjqKsrIyqFQq9s2dCEyqIt7BgMwxfD4fG3lHs+yGIh7RHXujF+9d8wBACYKbLEi3F+9c8d9of3oXqr92ecTnBqfOQo1xZ3reMrliZiyPSCSKOSoWhlLyy7wU3XhMDf59HEpLS1m32+TkJC/5YWZNiUxjNRqN0Ov1WL58edyRd6yiO/ZWL478x32C4KYCikb3Tb+E9bwOqx7YFv5pHN6P/j1v/cfy+PcwMBgMbJtF/0kQOTk5UCqV7DnS5UaLJaebKW0dgQwW3Wg9dbmKLkVRGBsbw/DwMJYsWTLLEstEqHxuysWCfyP23NxclJSUsLbFeBCLxZzvAsbfOYEjX7hfENwU0/eLv8FyTodPvHxvyPd5IkFAqB4GjFORiYqNRiOczovuOZVKxX6evF5vwvskXInlrtBqtQqim2649F+gKIoVM6ZDWajbLCZCTXRjIlbR9Z+yu2jRIqxfvx5erxdDQ0MJrSNalzCGifdO48jnfwjwZDQRiI2Jwz14o30nrnj/55AqAntCUxSVUJoqGJFIBLlcDrlcPmsShMPhwNTUFLxeL86ePQufz4esrKyA9IR/VMwXBEFwnvLCVHFkChkrutE20sKJrn/kWFRUhPb29ojf3slsZB4O/ym7/vPcCILgtctYOCbfP413//0e0ILgphXbgB6vN96AK9/7OVTLlrA/T1U1gVgsRm5uLrxeLwCgtrYWQGDzcZ1Ox242+/e7zcnJYfvdxgPXVBwTQMxVi3IoMlZ0gfDuKplMxt4eMdA0jYmJCWi12phG+KTSwmuxWNDf3x8wZdefZPXT9WfqX2fwzmfvAU0KgjsX8JkdONx2Gz75l/tQcvnFsqhUb2wFpzNCNR/3H+E+MzODkZEReL3egJE8TPNxLgIZax45k2qVM1p0w+EfnTKTgDUaDQoKCmIe4ZOKSDecZTeWY3AlUj/cqaNn8fZn9giCO8egfAT++cX7sfbxHai76d9TXjfLdT5aqBHu/lGxfwlmtKiYa6Qby4bbXGFeiq5MJoPX62XFNj8/H62trTFNAmZI9py0SJbdYJIZ6Ro6zwuCO5ehaZy863ewnB1G/u1Xp/R2OpHqhXBRscvlgs1mg9lsxujoKDweD6RSKSvCbreb0xdLvBbgdJLRohsqvUDTNCwWC6anpyGVSgNyovHAWC4Txf84XC27oY6RjJyuobsPb17536CJxNMoAslF+/wbUJ7oR8crP0zZOUmS5LVqIdzUYJ/Px0bFTqcT586dA03Ts6JiuVzOfl4yrUYXyOAm5sEwu/2dnZ0wmUxQqVRoaWlJeFeTr0hXKpXC6/Wir68Px48fR1FRETZu3BjTrDQ+bimDRdfQ3Yc3N98pCG4G4Tw1hH998i54Z2wpOV+q6nRlMhkKCgqwdOlSyOVytLW1oa2tDcuXL4dCoYDFYkFfXx+6u7tx/PhxPP/88/j973/P5pNj5dChQ2hoaEBdXR0ee+yxsM/r7u6GRCLBX/7yl0Quj2VeRLomkwmDg4NQKBRYvXo1FAoFPvzwQ17OIZVKZ/U7iBWCIDA1NQWTyYSGhgbOlt1k4C+6Q+8dx9F/vw8QUgoZh2fMhNdbbsTmtx9HflN8zV64kg4bMHDx88046FQqFZYs+biCw+fzwefz4cyZM9Bqtbjiiivgdruxb98+NDU1RT02l/lozPPuvvtufOYzn+HtujJadM1mM86fPw+5XB5yt58PEol0mf4No6OjKCwsRHl5OSoqKnheYWyIxWJ4vV7868XXMHzT04LgZjCkw423PrEbHc/fhYqtm5J2nnTYgKMhk8mwceNGTExMoLS0FA8++GBMqTcu89EA4IknnsCXvvQldHd387b2jE8vNDc3Y82aNUkRXCB+S/HIyAiOHTsGiqKwadMmlJWV8TbRIl7cbjf6+/thPDkI3befAQhBcDMdmqTw4Tf24tzjLyftHOmwAXPFP6crkUg4r5PLfDS9Xo9XXnllVqvHRMnoSLegoCCiICba1QuIbfOKMV5otdpZLje+6n3jwefzQaPRwGg0osgthemOfYCQw51XnPvxi7CcHcam3/8378dOtejGUhLHTNuOFS7z0Xbv3o29e/fyfu0ZLbqR4KOrF8At0g227La1tc2qBeZTdLl+mZAkCZ1Oh7GxMVRXV6PEK8Phz+4WBHeeov/7Ubx5yS5sfvcnkMj5qzZItejG2taxpaUl5nNwmY/W09ODa6+9FgBgMBhw4MABSKVSfPGLX4z5fP5ktOhGa3rj8/mSLrrhLLvB8D0RONKHgKZpjI2NYWhoCGVlZejo6IC9X4+Dl34XlJef8fQCcxPL2WEcaN6BK97/BZTlRdFfwIG5LrrxTALmMh9Nq9Wyf9++fTs+97nPJSy4QIaLbiSS7SSLZtnlepx41xPqQ+AfcRcWFrJ9JSx9Izi46XZQHn5G0wvMbTwGKw6tvRmXvfoQijdG38mPRqqnD8ci8vHW6XKZj5YsMlp0uUS6iRLclYurZTcYviPdYMxmM/r7+6FQKAIibsvAKA52CIK70KA8Phz5t3uw/onbsPz6KxM6Fk3TKRXdWBuYx9vWkct8NIbnn38+rnOEIqNFNxJ8RboMsVp2g+HaUjEaweJtt9sxMDAAiqLQ1NQU0CfVNjiGQxt2gnJ7Ez6vQAZC0ei97UlYzg5j7aPfSvdqOJMq0U0XguhGwe12w+Vy4fTp0zFZdpMFE+kyVmK73Y76+voAbzsA2DTjOLDhVpCC4C54Bn/9KiaP96H1j3chLy8voZaLqSAW0bXb7bOa7Mx1Mlp0k5le8J+yK5VKsWHDhjnRs1MkEmFoaAg2mw21tbUhp//ahyZwoO1WkC5BcAUuYjvWjw8v34PK330bhBRsy8Xc3NykNSKPF6453UzspQtkuOgC4XvqxmvfDTVlt6urK+WbCcFQFAWdToepqSlUVFSgo6Mj5Hocuim8vv4WkK7ErMsC8w/vqBG6r/0SV7z7E8irFrPNZYaHh+F0OkMOreRzQgVXCIKAUqnk/Py5HLWHIuNFNxyxRrqMZZdxqvhP2U10qKQ/sRo2/A0XpaWlqKioQFFRUWjBHZ3Ga+tuBukUBFcgNITViTc2fgefePlelF7ZGpCWCh7lrtFoQBAEPB4PNBoNGxVnZ2cnVei4ftZS3VeYLzJedCNFulxyuhRFQa/XQ6fTsTWtwbc2fJefcR1DYjAYMDg4iIKCArb8i6mcCMYxOo3X134bpNOd8DoF5jc0QeKDLz+E1Q9vx4rbP647DTXK3ePx4MyZM8jNzYXNZsP4+DjcbndA79vc3FzOEyG4wFV0bTZbwMZxppDxohuOaBOBI1l2g5FKpbzW2EZ7QzE1wHK5fJbhImQz9DEjXl93MwiHILgCHKGB0z94HpZzw2j/za6wT6MoCllZWVi8eDEWL17M/ty/9+3IyAicTifb+5aJiJmJELHCdRqE2WzOuF66wDwQ3XC3F+GGU3Kx7AbDZyPzSOLtcDgwMDAAkiTR2NgY8ls8uE7XNW7Ea2tvAmF3Jbw+gYXH8L53Ye0bxafffAziEEIXblOL6X3rX67F9LW12WwwGo0YHh6Gz+eDXC4P2LRTKBQR0wJc7wYzsYE5MA9ENxyhUgJcLbtcjhUP4UTX4/FArVbDarWivr4+YAx2pGO4xo14dc1NIGyC4ArEz0zvAA4078CVH/wS2YsDLbWxtHUMNSeNpml4PB42Kp6cnITL5YJEIgnYsMvJyWHPwzW9YLFY4rIAp5t5K7r+ZgR/y25LS0vMdX18To/wPw5BENBqtZienkZNTQ2ampqibgyIxWL4fD64Jmfw2tpvC4IrwAvuiRkcXHUjLj/8GArW1LA/T7SBuUgkQnZ2NrKzs1FcXMz+nCAIOBwONk9st9tBURSUSiVcLhdmZmaQm5sbsaZYiHTTRLTblBMnTsRs2Q2G70iX6bc7OjqKpUuXhi3/CncM97QZr239b/issY8oERAIB+ny4p3L70T7776Lqi9fevFnSWp2I5VKkZ+fHxCp0jQNp9OJkydPsgMrmTHu/nlipqY43mY36SbjRTcU/pbdlStXJmwTZOabJYpYLIbBYEBfXx9KS0uxcePGmMvQCIsDpz/3I5BChCuQBGiKQte3fgbreR1W/vA/U9phjKkTlslkqKurY3/u9Xphs9lgt9thNBrhdDrxq1/9CiaTCUuXLsX777+P1atXcxbgQ4cOYdeuXSBJEjt27MCePXsCHv/Tn/6EvXv3AgBycnLwm9/8BmvWrOHvOqP0A0i8WUCSoSiK3TALnrI7MDCADRs2JPymmZ6exszMDFasWBH3MQwGA86cOYOcnBysXr06rl1dt9GCf7TsAGFxxL0OAQGulH22Hct+th00TaOysjIl56RpGj09PWhvb4/4PJvNhoceegjAxaCor68Pb7zxRtT0HEmSWLFiRcBstBdffDFgTM/Ro0fR1NSEgoICHDx4EA888AA6OztjvZSwC5kXka6/Zbempoa1xg4PD3MuP4lEIukFJp+clZWFyspKyOXyuATXM2PDa6tuEgRXIGWMH+zGTJ8Oa/7y/ZSdk2sOOTc3FxKJBNdccw2uuOIKzsfnMhvtkksuYf/e0dGB0dHRGK4gOpllWg6B0+lEd3c3cnJysGnTJpSWlrLfduHKxmIlHtF1Op04deoUBgYG0NDQgDVr1kChUMRV7+u1OPDqqhtTNnJbQIDBrZlE76fvhUM3lZLzxTLtJZ6NNC6z0fx59tln8dnPfjamc0Qj4yNdpVIZYNn1J9mlXqHwer1Qq9WwWCyoq6sL2LGVSCRwu2MzMHgtDry68lvwGq0xvU4gdoo2NMDY1ZfuZcw5CIsDh9ffik/+9T6UXLY6uedKcltHLrPRGN599108++yz+OCDD2I6RzQyPtIViURhb0f4amTORbwJgsDg4CC6u7uxaNEibNy4MUBwgdgbmXttTry6cgc8BkFwU4GXFEYZhYPyEvjn1vuhfvZQUs+T7FE9XGajAcDp06exY8cO7N+/P2LdfDzMC9ENR7Lqa/1hun91dnZCJpOx49ZDrSsWO7HX5sSrLd+Cx2BJaO0C3PEZhfRNRGgaJ+54Gse/90zSThHLHozNZotZdP1no3m9Xrz00kvYunVrwHN0Oh2uueYa/OEPf0ho8zwcGZ9eABJvehONUFMfaJrG5OQkNBoNFi9ezKn8i6ud2Odw4bVVO+CZFgQ3lXhGDJCoskEKPSwiovmfg7CeG8anDvyY9y5fseR0Y3HLMXCZjfbQQw/BaDTi1ltvZV/T09MT24VEIONLxoCLedRQ1zE+Pg6Xy8XuVCbC0aNH2V1No9GIgYEB5OXloba2NmrvBgamt8LatWvDPsfncOHVlh1wT84kvGaB2Fm0ajnMH2mjP1EAyqWLceW/foms/MhDWWNhbGwMFEVFLVGjaRqXXXYZTpw4MVfbO4ZdVManF4DITW/4nJNmtVrR09ODkZERrFq1Cs3NzZwFF4ie0yVcHry26kZBcNOIjEcBme84R6ZxoOlbsJzX8XZMrjldJsiao4IbkXkhuuGI1t6RK06nEy6XC319faivr8fatWujjlwPRSTRJVwevLpyB1zjpkSXK5AAwky52CAcbrz1id3QH4jZPBD6eBxzuna7PSN76QLzXHQTrdP1er04f/48Tp06xfa2TcTrHU50CY8Xr626Ea4xY9zHFuAH+9BEupeQcdAkhQ+vexTnf/qXhI8VS1vHeHuppJt5Ibp8pxcIgoBarUZ3dzfy8/PR0dEBhUIRcmJDLITakGME16k3JHRsAX7wGqzIWpJ5navmAmcf/iOOffOnCR1jvrd1BOaJ6IYj1jpdpvtXZ2cnJBIJNm3ahPLycohEIt4amftzUXBvgnNkmtfjCiRGbvWSdC8hYxn92wd485PfBRXnHSZX0c3UqRHAPBHdcJFu8JSFcNA0jYmJCRw7dgxutxsbN27EsmXLAkwXfG/KkT4fXl9zE5wjqbFXCnBHnBV6bJMANywfafF60w64JmLfEOaa083UXrrAPBHdcHDZ2TSZTOjq6oLRaERrayvq6+tDftPyKbqkz4fX1nwbjmFBcOciHsFynTCeaQsOrLoR5197HxaLhbMpaL6P6gHmiTkiHmw2G/r7+yEWi7Fy5cqo1QixWnjDQflIvL72Zji0wobNXMWuGYdIJgHtS/z3vZChvQTOXv9zmO/7MrI+1QiSJNnBlUxT8uCOexRFceoyJohumokU0YpEooBfpMvlwuDgIFwuF1asWMH5F8dHpEuRJMZ3PAPfqFAWNpehPD7kNS6F9cJI9CcLRIamoX/w/6L+tq1Y89A2eDyeWYMrs7Oz2cGVFEWBpumod6lWqxXLli1LzTXwzLwQ3Uj4i6VGo8HMzAzb/SuWwupERZciSby+7mZBcDMEn3xeZ95SzsBT/4D53DA2/fleFBUVBXz+vF4v7HY7bDYbPB4Puru7IZPJWCHOzc2FUqkM+LxmcqQ7L95ZkcRTIpFAq9Wiu7sbubm56OjowOLFi2N2ssTSrCYYiiRxYP0tsA2E79spMLfIzuLuNBTgxvS7p/D2xl2g3V72DpSiKEilUixatAhLly6FUqlEW1sbWlpaUFRUBI/HA61Wi56eHvT09KCvrw9PP/00pqamOE/z9ufQoUNoaGhAXV0dHnvssVmP0zSN73znO6irq8Pq1atx/PhxPi49gHkb6VIUhbGxMczMzEChUKCjoyOhCRLxloxRJIkDbbfC2sdv93mB5OIWnIFJwaEZx8HmG3HV+79Abm05m05wu90YGBhAfn4+SJKESCRCbm4u8vLyIBKJWJFmIuK+vj7cfvvtkEql+PznP4/7778/6rlJksRtt90WMKpn69atAVMjDh48iIGBAQwMDKCzsxO33HJLPKN6IjLvIl2apjE1NYVjx47B6XSitLQUJSUlaRnZQ9M0Dm7YKeQGMxDXqAFZRZnpeJrr+KxOvHP1HngtDohEIkxPT+P06dMoLy9HU1MTsrKyIJPJIJFIWLElSRI0TUOlUmH37t1QKBR48803cfToUezYsYPTef1H9WRlZbGjevzZv38/tm3bBpFIhI6ODpjNZoyPj/N6/fMq0p2ZmUF/fz9UKhVaW1uRnZ0NtVqd9J66oaBpGgc23AbLueGEzy2QHnKWLYFJKB/jnUWrlmPD07tAZ0tx6tQpyGQytLW1QSa7WB/NbHr7B0pMvb3b7cbPf/5zjIyMsPMGKyoqOJ031Kie4Cg23DifsrKy+C42BPNCdGmaxqlTp0BRFFpaWpCTk8M+lo45aTRN47X2W2E7KwhuJiNRCHldPhHJpCi5YTMUX27HGesYiKM6LFmyBKWlpVFfKxaLcfLkSezatQtbt26FVqtlRZorXEb1xDLOJ17mheiKRCI0NjaGbLMok8ng9SbeOYrrRprT6cShS74Dd99YwucUSC8+qzB5mS8K2xuw4dffQdayxTh//jxy8/NQXl4Op9OJ6elp9o7Uv443NzcX2dnZ8Hg8ePzxx3HkyBE8++yzWL06vjltXEb1cB3nkwjzQnQBhG1II5VK4XAk/uGJtpHGNMn56KuPw9vPbw5IID3YNROASAREbvQvEAGJUo5V912P+ls+D/3YGPSnTmHFihUoLCwEgICyL5qm4XK5YLVaYTabcfjwYTz66KPwer1YsWIFdu/ejSVL4u+L4T+qp6KiAi+99BL27dsX8JytW7fiySefxLXXXovOzk7k5+fzmloA5pHohoPPkT2hRJ2iKIyOjkKn08F0x58EwZ1HEHYXcmrKYNcIv9N4KLl0Fdqfuh3i0nycOHkSOTk5aG9vD7upLRKJoFQqoVQqUVBQgH379qGsrAwPPvgg3G43jh8/jqKioriFl8uoni1btuDAgQOoq6uDUqnEc889F/f1h2NejOsBLkaaoW7/bTYbtFpt3Lck/viP7KFpGtPT0xgcHERxcTE0O34NU3d/wucQmFsUX9IMw9Fz6V5GRiHNVWDNw99EzQ2fwcjICMbHx9HY2MjZzNDb24vdu3fjK1/5Cu68807OM9PmGGETwRl5NbHAd3cw4KIbpq+vD0qlEmvXrsX7W34gCO58JWPCjrlB2dXr0far20AtUqC3txcFBQURo1t/3G43Hn30UXz44Yd44YUX0NLSkoIVp555L7p8jewBLqYSTp06BZ/Ph8bGRqhUKrx91d0wdl7g5fgCcw/3hGCS4EJWYS7WPrYD1ddejuHhYUyf06KxsZHzdIfu7m7ccccd+NrXvoYjR45kanTLiXlzZeHKOvhoPu7z+aDRaOByuVBfX4+SkhKQJIm3P7MHhg+FW8/5jH1oEpJcBUibK91LmbNUbt2E9b+4BT7lxVHlxcXFaGtr49QtzOVy4ZFHHkF3dzf++Mc/oqmpKQUrTi/zwpEWiURq7CiKwvDwMLq6uqBSqZCfn4+8vDwQBIF3t9wLw7/O8rhSgTkJTUNUkpkDEJNNdskiXPKHPdj0xz0YtRlx4cIFNDc3o6amhpPgdnZ24uqrr0Z5eTneeeedBSG4wAKIdOOBsRKr1WqUlJRgw4YNEIvF8Hg86O3txfSeF+Hq1fJ2PoG5jap4ESxqoeG8PwVbWtH44DcgK8hBV1cXSktLsX79es7R7Y9+9CMcP34c+/btQ0NDQwpWPHeY95EuQ5QqDRaLxYLu7m5MT0+jtbUVy5cvB3Ax6q2urobn0QOC4C4wJD5hN41BUVGMTS/fg3VP7cSkfQYXLlwARVGYnp5Gf38/xsbGYLfbw47J+vDDD3HVVVehqqoK77zzzoITXGCBRLqMmyxSct7pdGJgYAAEQaCxsRE5OTls6zmmy9F7X7wfk++cTMbyBeYwwlh2ACIRarZfjbU/vgF20oO+vj5UVFSgsrISIpEIPp8PNpuNLdF0OBwQi8XIzc3FmTNnUFxcjAMHDuDs2bN46aWXsGLFinRfUdqYN3W6NE2HtfseP34cTU1NIftvMptkJpMJ9fX1KCoqYjsaMWILAP/88oMYP9ST1GsQmLsoly5esFObVctL0f7k7Sj6RDMGBwfhdDrDfp78IQgCNpsNv/jFL/DGG2/AbDajrKwM69atw69//WtOqYgMZuHW6QKha3WZceujo6Oorq5GfX09aJpmbmYufQAAHk5JREFUn+cvuO9/9WFBcBc4ivLiBSe6IrEY9bd8Dqvu/wYsLge6u7tRVVWFhoYGTnsoHo8Hjz76KM6dO4e//e1vqKurg9PpRF9f33wX3IjMmyuP9Cbwr9WlaRqTk5M4duwYfD4fNmzYgLKyMpAkyaYSxGLxx4J77Y8wdqArJdcgMHcRSefNR4UTeQ1LsfnNvVj5o+3oH9JAp9Nh3bp1qKioiCq4NE3j/fffx1VXXYXGxka8+eabqKurAwAolUqsW7eOlzXecMMNKCkpwcqVK8OuI9lTIOJhXkW6IpEo5IYZ097RbDajv78fSqUSra2tkMlks/K2/vzrPx/B2Gv8do0XyEy8hoXRV1cklaBx9zVo+f51MFnN6OnpwbJly1BaWsopurXb7bj//vvR39+Pv/71r6itrU3aWrdv346dO3di27ZtIR9PxRSIeJhXohsOmqahVquRlZWFpqYmqFSqiGILAEe/8RhG93+YhtUKzEVsmnGI5TJQnsR7M89VFq2uwYbffAeqpqU413cBNE2jtbU1ZMvUYJjods+ePbjpppvw1FNPJT2FcNlll2FoaCjs4+GmQPDdNSxW5pXoBke6Pp8ParUak5OTWLx4MZqamkCSJDuDKdyb4uh/PY6RV/6VqmULZAC0j0BuSzWs87AxvVguQ8vdX0PjHV/GtNGAC729qKmp4dzNy2az4b777oNGo8Hf//73OTMaPRVTIOJhXiaqKIrC0NAQurq6kJuby5aneL3eWVUJwZz4/v9g5G/vp3K5AhlC1qKc6E/KMArbG3D1B79E7a4v4qOzZzA9PY22tjZOgkvTNI4cOYKrr74a69atw+HDh+eM4AKpmQIRD/Mq0qVpGhMTE9BoNFiyZAk2btwIkUgEt9uN6elp9Pb2QiQSIS8vD3l5ecjPz4dSqWR/ETRNo+RTqyGSy6D/4DSc50dBWZxpviqBuQLt47dbXToRZctQ9K3NqPzmlRiHE8aeHravCBdsNht+8IMfQKfT4R//+Aeqq6uTvOLYScUUiHiYV6I7ODgIt9sdsElG0zTkcjnbJo4kSVitVlgsFqjVajgcDmRlZSEvLw85OTmwLcuDe1EDNty2BcXFxXBoJ2Ds6Yfp+ABMvf2YOakG6Up8/I9A5uGYJyVjTHNxqliF8+fPA7hYVaBWqzEyMsKOPs/Ly4NCoZg1bfvIkSO45557sHPnTjzzzDNztvwrFVMg4mHemCOAi+kDpvQrWhrBH4/HA61Wi4mJCWRlZbEd7JloOC8vj3WzUQQJy7lhmHr6Yezth6lnANYLOtBkaNujwPwiu2QR3FPmdC8jLmR5Sqx+eDtqvvkZjI2NYXR0lDUEMXi9XlitVlitVthsNjidTmRlZaG7uxsURaG7uxsWiwW//e1vUVVVlb6LAXDdddfhyJEjMBgMWLJkCR588EF2CO3NN98Mmqaxc+dOHDp0iJ0C0dbWlqrlhRWeeSW6d911F3JyctDW1ob169cjNzc3quiaTCYMDg6ioKAAy5Ytg0wmY2c1WSwW9g1IkiRycnJYEc7NzWW/4QmHG6YTgzD19sP0/6Nix7DQIGU+ktVYBu+FzBvfwzQXFxXl4Ny5c8jJyUFdXR2n5uIejwe/+93v8Morr8Dj8YCiKBQWFuLPf/4ziouLU7D6jGRhiG5fXx+OHTuGzs5OHD9+HF6vFytXrsT69evR3t6OlpYWdmzzzMwMhoaGIJFIUF9fH9XSSFEU7HY7m5qw2+1sfpgRYv/8sHvawoqwsbcfpt4BeE22pP8bCCSXxZ9owXQGtfRkmosvu+7TGBkZwdjYWEyjcywWC+655x5MTU3h6aefZqsBpqenUVhYyEm0FygLQ3SDcbvdOHnyJI4dO4bu7m6cPXsWMpkMMpkMcrkcP/nJT9DY2Bh3ToogCDYStlgscDqdkMvlAULsX+No14xfzA/39sPYOwDzKSE/nGnkr1wGy5mhdC+DE5VfuATrf34zyJwsnD9/Hvn5+aipqeEklDRN480338R9992HO+64A9u2bUta7vbQoUPYtWsXSJLEjh07sGfPnoDHLRYLrr/+euh0OhAEgTvvvBPf/OY3k7IWHlmYohvMX//6VzzwwAPYsmULsrOz0dPTg+HhYVRWVqK9vR3r169HW1sbCgoK4i4t8Xg8bFrCYrHA6/VCqVSyIhycHzaeVqP/0FFYTmpAaw1w9OuF/PAcRqKUg/T4gDn8O8ouWYTWn92Mii9sgk6nw+TkJBobG5Gfn8/p9WazGd///vdhMpnw9NNPo6KiImlrJUkSK1aswJtvvsl+Dl988UU0Nzezz3nkkUdgsViwd+9eTE9Po6Ghgd1/mcMs7IY3DB0dHejq6gpIJTA1vZ2dnThy5Ah++tOfwmazoampiRXhNWvWIDs7m9M55HI5SkpK2NIbmqbhdDphtVoxNTWFwcFBUBSF3Nzci5UUTitqtl+NsrIyiEQiIT88xyGdHuTUV8A+oE/3UkJSfd2nsW7vjfDKgJ6eHhQVFaG9vZ1TlErTNA4fPowHHngAd955J66//vqkVyZ0dXWhrq4ONTU1AIBrr70W+/fvDxBdkUgEm80GmqZht9tRWFiY0TPUMnflcRDqG1ssFqOmpgY1NTW47rrrAFx0sn300Ufo7OzE73//e5w+fRpSqRStra1obW1FW1sb6uvrOd2miUQiqFQqqFQqtlzFYDCgv78f2dnZyM/Ph06ng16v/7h+uLUGiz/RIuSH5yjZi/PnnOgqKorR9qvbUHpVK7RaLYxGI5qampCby23U0MzMDPbs2QOr1YpDhw6lrJ41lGssuD/Czp07sXXrVpSXl8Nms+Hll1+es2VqXFhQossVmUzGCuwtt9wCmqZhs9nQ09ODzs5OPPzwwxgcHERJSQkbDTMunkhpCZfLhf7+i6Pa165dC6VSyT7GNIG2WCyYmppi88NMWqJ48xqU/1s7+3whP5w+yLlkkvBrLu4Cge7ubpSUlHAeDEnTNA4ePIgHH3wQd999N77+9a+nVNC4uMYOHz6MtWvX4p133oFarcZVV12FSy+9lPOk4bmGILocYKoUNm/ejM2bNwO4+GYZGxtDZ2cnOjs78Zvf/AYGgwH19fVsyVprayuUSiWsViuGh4fZacKFhYWzziGTyVBYWBjwmNvthtVqhdlshk6ng9fr/XhAZkEeKr/0SVR/9VMAhPrhVGIbmRvpHmlZAarv/ypKL1+L/pGL0xpWrlwJlUrF6fUmkwl33303XC4X3njjjbQYB7i4xp577jns2bMHIpEIdXV1WL58OS5cuIANGzakerm8sKA20pINSZI4f/48Ojs70d3djd7eXhgMBhAEgW3btuHzn/88mpub485H0TQNh8PBbtIxeS7GQZSfnw+VSiXUDycbkQiyPCV8Fkd6Tv//m4s33ft1TFtM0Gq1yMrKAk3TkEqlbJoquIyRgaZpvP7663j44Ydxzz334Nprr01bTwKCILBixQq8/fbbqKioQHt7O/bt28c6SAHglltuwZIlS/DAAw9gcnISra2tOHXq1FyvERaqF1INTdP43Oc+h+rqamzZsgV9fX3o7OzEhQsXkJ+fz9YOt7W1oaKiIu5bOpIkYbfb2YoJu90OiUTCpiXy8/ORnZ3NfqhMw+M4+9p7IAYmINKaMHNiUMgPx0HB2lrMnFSn/Lx5DUvR/uvvoKCtHoODg7Db7WhqamJTVUyaiilldDgckMlkyMvLQ1dXF2pra/Hb3/4WJEniqaee4txJLJkcOHAAu3fvBkmSuOGGG3Dvvffi6aefBnDRWTY2Nobt27djfHwcNE1jz549uP7669O86qgIopsOLBbLrDIdmqZhMBjYtERXVxf0ej2WLVvG5oZbW1uRn58fd/Th8/nYaNhqtcLlckEul4MkSXi9XjQ0NARECUJ+OHYWf3Ilpj84k7Lz+TcXtzhs6O/vR2VlJadJDl6vFxaLBXv27EFXVxcbXW7evBnf//73U3QFCw5BdOcyFEVhcHCQFeHe3l44nU60tLSwQrxy5UpOzaSD8e+8VlBQAKlUCqvVCoIgoFKp2Gg4NzeXrcagCBKWs0Mw9Q4I+eEwFKyrxcyJ1ES6THPx3JZq9Pf3w+PxoKmpiXMZo8FgwPe+9z2IRCI8+eSTKCkpwcTEBIaGhtDR0ZHk1S9YBNHNNLxeL06ePMkK8ZkzZ5CdnY1169axQlxTUxMxLWG329HX1welUona2tqAYnImP8xEwzbbxRRDcH6YiaKC88PG3gE4dQs3PyzLVyU9p+vfXNxknsHAwACqq6vZmu5o0DSNv//973jsscfwwx/+EF/5yleSmruN5iwDgCNHjmD37t3w+XwoLi7Ge++9l7T1pBlBdDMdmqZhNpvR3d3NbtRpNBpUVFSgtbWVddQVFxfDZDJBrVZDJBKhoaGBc2kNSZJs2RqTD5RKpQH5Yblc/nH98JT5YsvLBVo/rKwugTNJG5NFGxqw4de7oKgtRV9fH0iSRGNjI+e7nampKXzve9+DTCbDE088gcWLFydlnQxcnGVmsxmXXHIJDh06hKqqKkxNTXHu35uBCKI7H2HGyB87dgxdXV3o6urC0NAQaJrGV7/6VWzZsgVr166d1RM1FphWf4wQu91uKBSKgP4STBMhYGHlh+UrK+E5M8rrMSVKOVbddz1W3LoV0wYD1Go1ampq/l975x4UZdn+8e+zHF6w5XxICUVkBXaVReQwxCQvL4wSvGaGls34k0wpTQmaGmb62fSLcUxEwsZXzEMhpolYVuPrG8ivUTJSOW0myiKsooHAqhwFhGVx7/cP2sfdOD0Cu8tu92eGmT3c7H0jD5f3872/13XB1dWV8+72u+++w86dO5Gamoq4uDi9OBMuXbqE1NRUFBYWAgDS0tIAQEsz/uyzz9DU1IRt27bpfD1TAJoGbIrweDx4eHjAw8MDq1atwqpVq+Dv749Vq1ZBKpUiLy+P9TcuWLCATeTw8fHhXB3K0tISzs7O7MEbIQR9fX3o7OxEa2sr6urq2LKXtra2sHOyw8yVi1j/sFLRD+n/X0BLeQ2smnvQVXnbZPRhvo0tFJP4eeri4pbuTrh67RoYhkFgYCDnGgN3797Fe++9B2traxQVFenVUsUls6y2thZKpRIRERHo6upCcnLyiJ18TRkadE2I7Oxs8PmDfbwWLVqEDRs2sPnqEokEpaWl2LFjB2pra+Hk5ITAwEAEBgYiJCSEc4tthmFgbW0Na2trTJ8+HcDgjlutDzc2NqKrqwsMw8DCwgLd3d2YIZwJv3/+fUT/sLHqw/0tnZPyOZrFxe/evYvbEgkEAgFnSUClUuHbb7/FJ598gq1bt2L58uV6991yySwbGBiARCLB2bNn0dvbi2effRahoaFsD8O/CkYddMcS7gkhSE5ORn5+PqZNm4bDhw9j4cKFBlqt7lEHXE0YhoGNjQ0iIiIQEREB4LGjoaysDCUlJcjOzoZcLodAIGD9wwEBAeDz+Zz+eHk8HmxsbNg8/97eXlaHdHNzQ09PD0pLS9m2SHZ2drANEsD1ufnsZxijPtxd1wyetSVUE5BPnCP9EfSvzfjbdAdUVlbCwsICQUFBWpLNaMjlcrz77ruwsbHBTz/9pNUFQp9wySxzd3eHs7MzW4skPDwcV65c+csFXaPVdLkI9/n5+dizZw/y8/NRWlqK5OTkIbc8lEEePXqE2tpaVh++fPky+vv74efnxwZikUg0ajBQqVSor6+HXC4f0gYGeOwXVWvECoVixLZIgHHow+Otr2vhwMecLS/D6h8itLS0QKFQwN7eHi4uLmxnktEkIJVKha+//hqffvoptm3bhmXLlhm00y2XzLLq6mokJiaisLAQ/f39CAkJQV5eHubPnz/KJxstpqfpcikJd+rUKcTHx4NhGISGhqKjowPNzc1TojndVMPMzAxCoRBCoZAtEN3X14fLly+jpKQEe/bsQVVVFfh8vlaRn1mzZoHH40Emk6GtrQ0uLi4ICQkZ1spmaWkJFxcX9rZZsy3S/fv3cfPmTe22SE62Wvqw2j/MNgqdAv7h8bRlVxcXJ7ZWqK6uhrOzM+bMmcPWYm5qahrRwsfj8SCXy5GcnAxHR0ecP39+2Foe+sbc3BxZWVmIjo5mM8vmzZunlVkmFArx/PPPQywWg8fjISEhwVQD7qgY7U735MmTOHPmDL744gsAwNGjR1FaWoqsrCx2zNKlS/H+++/jueeeAwBERUUhPT1dn83pTApCCFpbW1FeXs5246irq2MljC1btky4CPyf2yJ1dXWBx+ON2BbJ0PqwY5A32ipqOY21crXHwl0b4f5iGO7cuYPGxkb4+PjAwcFh2PFqC5/632Lnzp2QSqVoa2vDmjVr8MYbb2Du3LlGXebQhDG9nS4X4Z7LGAp3GIaBs7MzYmJiEBMTA4lEgvXr12P9+vVwdHREUVERMjIy0N3dDZFIxO6IxWIx5+wpdYC1tbWFu7s7AO22SDdu3GA71NrZ2Q1+Bc/V0oe7m1sg/U8xOi7fgPmdTnT+VqczfbiHY4BXFxcf+BsPEokEtra2CA4OHlVCMDMzg729Pezt7dHc3AyFQoGQkBC89NJLqKmpwZYtW7Br1y54eHhM1o9D0QNGG3S5CvdjjXkSxjq4O3bsGNLT0wEMHmrt27cP/v7+455vqiMSiVBcXMweoK1evRrAoHarLgKfk5ODq1evwsLCAgEBAaw+LBAIOO/QzM3Nh5S9VCgU7A6woaGB1YfNzc3R3t4Oj+eDEJzw+BRfV/qw4l4HrKY7oE/ePuz76uLiM5YE4vfff3/i1jkqlQq5ubnIyspCWloaYmNjdbZx4JJRBgDl5eUIDQ3FiRMnsHLlSp2sxZQxWnmBi3D/ww8/ICsriz1IS0pKQllZ2bjm43Jwd/HiRQiFQjg4OKCgoACpqan04A6DdxwPHjxgi8CXlZXh5s2bePrpp7X0Ya4JAMPR19cHqVQKpVIJPp+Pnp4eti2SWpbg8/lsoNfShyUytEnGrw87BfugtbxG+0WN4uIKngrV1dVwcHAYM3Vbk6amJiQlJWHGjBnIzMzk3MF3PHC5vtXjFi9eDCsrK6xbt44G3ZExPXmBi3AfGxuL/Px8CAQCTJs2DTk5OeOej8vBXVhYGPs4NDQUd+5MbraSscIwDOzs7BAVFYWoqCgAg4G4sbERpaWlKCkpwd69e9Ha2gpvb2+2CHxAQMCw9WA1UX9OQ0PDEG+rSqViNdH6+nq27CXbFmnuDHiJ50CwPgbA+PVhnqX2n9FTntMRnPU2XBbNx+3bt9HS0vJErXNUKhW++uor7Nu3D+np6YiOjta5LMbl+gaAPXv2YMWKFSgvL9fpekwZow26ABAbG4vY2Fit1zZu3Mg+ZhgGe/funZS5uGTcaJKdnY2YmJhJmdsUYRgG7u7ucHd3x4oVKwAM7qKkUilKS0tx8uRJfPDBByCEQCwWs7thX19f1lbW3t6OmzdvwsbGBsHBwUOKw/N4PFb3Vf/uxmqLZB/iPdQ/LKl97JgYxj/c39E9+DPxeJi76QX4/d//4OFAPyoqKuDi4sK5dQ4wKIElJSVh5syZ+PnnnznLEBOFy/Xd2NiI77//HufOnaNBdwIYddDVJ09yKFdUVITs7Gz88ssvul6WSWFmZgY/Pz/4+fkhISGB7aQskUhQVlaGjIwM1NTUsIV3Hj58iH379kEgEHDeCT5xWyRbWzy9JBBuMY9bw/xZH35wvQF282cj6F+b4Rjkjbq6OrS3t0MkEg2bsDIcKpUKR44cwYEDB5CRkYHFixfr9dCXy/X9zjvvID09nXMKOWV4aNDlCNdDucrKSiQkJKCgoMBg2UGmgrqTcnh4OMLDwwEA165dw5o1a+Dl5QU3NzekpKSgqakJnp6eWkXgbW1tOQctKysrWFlZsRWv1MG+s7MTcrkcMplMuy2Si90Q/zBA8KC7G2VlZZgxYwaCgoI4z9/Q0IC3334bc+bMQXFxsUEaLnK5visqKvDqq68CGKzRm5+fD3NzcyxfvlyvazV2jPYgTd9wObirr69HZGQkjhw5oqXvUiaPe/fuoaurC15eXuxr6iLw6mw6iUSCvr6+IUXguRaOGY7R2iLx+Xy0t7ejt7cXIpFIq8vzaKhUKhw+fBiff/45MjMzERUVNaV7lWmydu1aLF26lB6kjYzpHaTpGy4Hd1u3bkVrays2bdrEfk9FRcW456QWnqG4uroOqcHK4/Hg7e0Nb29vtmqVQqFgi8AfOHCALQK/cOFCNhB7enpy1lrVAVZTY1UqlWhsbIRMJmMDek1NjVYix0iBvr6+HomJifD29saFCxc4yxC6gsv1TZkc6E53ikItPJOLugh8WVkZWwT+1q1beOaZZ9ggHBgYCCcnJ067zYGBAchkMvT29kIoFMLa2hqEEDaVV+0h1myLpPbonjhxAjk5OcjMzERkZCRN2DFN6E7X2KAWnsmFYRg4ODggOjoa0dHRAB4X6CkpKUFxcTF27dqFzs5O+Pr6skkc/v7+sLa21vqslpYWtnWOr68vGzQZhmH1YXWXXc22SPv378elS5fQ19eHF154AfX19VAqlROSPSjGBw26UxRq4dE9PB4Ps2fPxuzZs9kDIqVSiaqqKpSUlODYsWNISUkBj8dDQEAAfH198eOPPyI+Ph7R0dGcUpvV9Ydzc3Nx/fp1fPnllwgODsaVK1dQUVExxOZGMX1opYwpCrXwGAYLCwssWLAAGzduRHZ2NioqKnDu3Dl4eHggIyMDKpUKH3/8MVauXImPPvoIp0+fhlwuH/b3BQC3bt3Ciy++iFu3buHChQuIiIjAU089hbCwMCQlJU16sZozZ87Ax8cHAoEAO3bsGPL+sWPHIBaLIRaLERYWhitXrkzq/JSxof/NTlGohWdqwDAMmz5cWVkJFxcXEELQ3NzMFoE/ePAg7t27xxaBDwoKgr+/P44fP46jR49i9+7dWLRokc6120ePHmHz5s1a5wDLli3TkqQ8PT1x/vx5NlX9zTffpKnq+oYQMtoXxUAolUri6elJ6urqiEKhIGKxmFy7dm3E8a+99hr55ptvJjxvQUEB8fb2Jl5eXiQtLW3YMUVFRcTf35+IRCISHh4+4TlNgYGBAVJVVUUOHTpENmzYQDw8PMjLL79Menp69LaGixcvkiVLlrDPt2/fTrZv3z7i+La2NuLm5qaPpf0VGTGu0p3uFMUQFh4uO6WOjg5s2rRJq402ZdBSJhKJIBKJ8Prrr4MQondXAk1VNw5o0J3CjFVbQpPDhw9PeD4ujonc3FzExcVh1qxZADDEM0sZxBA2MEJT1Y0CepBGYRlup9TY2Kg1pra2Fu3t7YiIiEBgYCCOHDmi72VSRuBJU9VPnTpFU9UNAN3pUli47JRoG+2pS3BwMGQyGZv0kZeXh9zcXK0x9fX1iIuLw9GjR+nvzEDQoEthoW20jRtDpKpTxsFop2z6P/CjGBIujgmpVEoiIyOJUqkkPT09ZN68eeTq1avjnnMst0RHRwdZunQpEYvFRCQSkUOHDo17LgpFj1D3AmVs9N1Gm4tbYu/evRCJRDh9+jTu378PHx8frF69mqbOUowWGnQpWnBxTKSkpCAlJWXCc3FxSzAMg66uLhBC0N3dDUdHR5o6SzFqqHuBYjC4uCUSExNRXV0NNzc3+Pn5Yffu3ZOeOjuVGCuNlxCCpKQkCAQCiMVi/PrrrwZYJWUimO7VS5nyEA5uicLCQixYsABNTU347bffkJiYiAcPHuhriXpFLbcUFBRAKpXi+PHjkEqlWmMKCgogk8kgk8lw8OBBvPXWWwZaLWW80KBLMRhc3BI5OTmIi4sDwzAQCATw9PTE9evX9b1UvaApt1haWrJyiyanTp1CfHw8GIZBaGgoOjo60NzcbKAVU8YDDboUg6HpK+3v70deXh6WLVumNWbWrFk4e/YsAODu3buoqalhNeDxsm7dOri6uo54AGioW3gucguXMZSpDQ26FIOh6ZYQCoV45ZVXWLeE2jHx4Ycf4uLFi/Dz80NUVBTS09Ph7Ow8oXnXrl2LM2fOjPi+oW7hucgtXMZQpjZjteuhUEwShmFmA/gPIWTIdpdhmAMAfiKEHP/jeQ2ACEKITu/jGYZ5FkAqIST6j+f/CwCEkDRDr40yedCdLoUylGcANGg8v/PHa7qmHMBchmE8GYaxBPAqgH//acy/AcQzg4QC6KQB17ighkcKZSjD3a/r/JaQEDLAMEwigEIAZgAOEUKqGIbZ+Mf7+wHkA4gFcAPAQwCv63pdlMmFBl0KZSh3AMzUeO4OoEkfExNC8jEYWDVf26/xmADYrI+1UHQDlRcolKHQW3iKzqA7XcpfDoZhjgOIAODMMMwdAB8BsADoLTxF91D3AoVCoeiR/wItJizHnFyNtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "ax.plot_trisurf(x, y, z, triangles=tri.simplices, cmap=plt.cm.Spectral)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87ac237b",
   "metadata": {},
   "source": [
    "### Remeshing\n",
    "!pip install meshio <br>\n",
    "https://www.cgal.org/ <br>\n",
    "https://www.meshlab.net/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2d1c8a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygalmesh\n",
    "import meshio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "e177dc05",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = pygalmesh.remesh_surface(\n",
    "    \"bunny.obj\",\n",
    "    max_edge_size_at_feature_edges=0.0025,\n",
    "    min_facet_angle=25,\n",
    "    max_radius_surface_delaunay_ball=0.1,\n",
    "    max_facet_distance=0.001,\n",
    "    verbose=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "d71506da",
   "metadata": {},
   "outputs": [],
   "source": [
    "meshio.write('bunny_remesh.obj', mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "34848ddc",
   "metadata": {},
   "outputs": [],
   "source": [
    "v,f = igl.read_triangle_mesh('./bunny.obj')\n",
    "v2,f2 = igl.read_triangle_mesh('./bunny_remesh.obj')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "402f170a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8b4ff3066f934c6b9d5999c9f66c1a3b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016840…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "df06a10fee764c3088567adf9bbedf6d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016726…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb181c46f98>"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(v,f)\n",
    "plot(v2, f2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "fce35467",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total vertices :  34834\n",
      "Irregular vertices  no :  3147 9.03%\n",
      "Area average: 5.697847766042785e-06 Min:1.8044666802109844e-12 Max:0.00012397079397837317 Sigma:2.261048632168552\n"
     ]
    }
   ],
   "source": [
    "irregular = igl.is_irregular_vertex(v2,f2)\n",
    "v_count = v.shape[0]\n",
    "irregular_v_count = np.sum(irregular)\n",
    "irregular_ratio = irregular_v_count / v_count\n",
    "print('Total vertices : ', v_count)\n",
    "print('Irregular vertices  no :  {} {:.2f}%'.format(irregular_v_count, (irregular_ratio * 100.)))\n",
    "\n",
    "area = igl.doublearea(v2,f2) /2.0 # igl.doublearea computes twice the area for each input triangle\n",
    "area_avg = np.mean(area)\n",
    "area_min = np.min(area)\n",
    "area_max = np.max(area)\n",
    "area_ns = (area - area_avg) / area_avg\n",
    "area_sigma = np.sqrt(np.mean(np.square(area_ns)))\n",
    "print('Area average: {} Min:{} Max:{} Sigma:{}'.format(area_avg, area_min, area_max, area_sigma))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "283469e0",
   "metadata": {},
   "source": [
    "## Tetgen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3be4a53d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista as pv\n",
    "import tetgen\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "5a3a2482",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object `sphere.verts` not found.\n"
     ]
    }
   ],
   "source": [
    "sphere.verts?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "81cb0bf8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f2ff82369a6f414f845bf2d3c11a5438",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sphere = pv.Sphere()\n",
    "tet = tetgen.TetGen(sphere)\n",
    "vv, tt = tet.tetrahedralize(order=1, mindihedral=20, minratio=1.5)\n",
    "grid = tet.grid\n",
    "grid.plot(show_edges=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "491f0cd3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3298, 3), (17634, 4))"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vv.shape, tt.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "531b46e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(88170,)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid.cells.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "86a61552",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(17634, 4)\n"
     ]
    }
   ],
   "source": [
    "cells = grid.cells.reshape(-1, 5)[:, 1:]\n",
    "print(cells.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "a20eed0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2240, 3)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sphere_mesh = tet.mesh\n",
    "faces_sphere = sphere_mesh.faces.reshape(-1, 3)\n",
    "faces_sphere.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "0c657f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "verts_sphere = sphere_mesh.points\n",
    "# faces_sphere = sphere_mesh.faces "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9c8a26b9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((842, 3), (2240, 3))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "verts_sphere.shape, faces_sphere.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "cf7163b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "v = verts_sphere.copy()\n",
    "f = faces_sphere.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "709fcd5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import igl\n",
    "from meshplot import plot, subplot, interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "ad445842",
   "metadata": {},
   "outputs": [],
   "source": [
    "v,t,f = igl.read_mesh('bunny.mesh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "ec17a8e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bbfb0a05740d47ff99e1f8f850507824",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016860…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb184660b38>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(v,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "368eecbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "v1 = 0\n",
    "v2 = 1\n",
    "v3 = 2\n",
    "v4 = 3\n",
    "v5 = 4\n",
    "v6 = 5\n",
    "v7 = 6\n",
    "v8 = 7\n",
    "\n",
    "bbox_tri = np.array([[v1, v2, v5], [v2, v5, v6], [v2, v3, v6], \\\n",
    "                [v6, v7, v3], [v3, v7, v8], [v8, v4, v3], [v4, v8, v5],\\\n",
    "               [v5, v4, v1], [v5, v6, v8], [v8, v6, v7], [v2, v4, v1], [v2, v4, v3]\n",
    "               ])\n",
    "\n",
    "zmin = v[:,-1].min() + 0.1 * v[:,-1].min()\n",
    "zmax = v[:,-1].max() + 0.1 * v[:,-1].max()\n",
    "ymin = v[:,-2].min() + 0.1 * v[:,-2].min() \n",
    "ymax = v[:, -2].max() + 0.1 * v[:, -2].max() \n",
    "xmin = v[:, -3].min() + 0.1 * v[:, -3].min() \n",
    "xmax = v[:, -3].max() + v[:, -3].max()\n",
    "\n",
    "bbox_verts = np.array([[xmin, ymin, zmax],[xmax,ymin,zmax],[xmax, ymin, zmin],\\\n",
    "              [xmin, ymin, zmin], [xmin, ymax, zmax], [xmax, ymax, zmax],\\\n",
    "              [xmax, ymax, zmin], [xmin, ymax, zmin]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "e2e85d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_new = np.concatenate([v, bbox_verts], 0)\n",
    "f_new = np.concatenate([f, bbox_tri+len(v)],0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "6b687774",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ea9ec6a2c903446e88bee7470157540c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(0.0089205…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb1fcb69278>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(v_new, f_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "9eada13b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tetgen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "10ba8571",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9617cb7c9c634b5384e6ce0e0f3cb34c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.016860…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f560e4c75d57426b98cb4a4fab6565d8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "ViewInteractiveWidget(height=768, layout=Layout(height='auto', width='100%'), width=1024)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v,t,f = igl.read_mesh('bunny.mesh')\n",
    "plot(v,f)\n",
    "tet = tetgen.TetGen(v, f)\n",
    "tet.make_manifold()\n",
    "vv, tt = tet.tetrahedralize( mindihedral=2, minratio=1.5, verbose=2)\n",
    "grid = tet.grid\n",
    "grid.plot(show_edges=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "66295b72",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista as pv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
    "\n",
    "from skimage import measure\n",
    "from skimage.draw import ellipsoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "12ec3a67",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "# Generate a level set about zero of two identical ellipsoids in 3D\n",
    "ellip_base = ellipsoid(6, 10, 16, levelset=True)\n",
    "ellip_double = np.concatenate((ellip_base[:-1, ...],\n",
    "                               ellip_base[2:, ...]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "29186967",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(15, 23, 35)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ellip_base.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "c2570b0b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(27, 23, 35)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ellip_double.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "e0e5996b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: seaborn in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (0.11.1)\n",
      "Requirement already satisfied: pandas>=0.23 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from seaborn) (1.1.5)\n",
      "Requirement already satisfied: numpy>=1.15 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from seaborn) (1.19.4)\n",
      "Requirement already satisfied: scipy>=1.0 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from seaborn) (1.5.3)\n",
      "Requirement already satisfied: matplotlib>=2.2 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from seaborn) (3.3.4)\n",
      "Requirement already satisfied: pillow>=6.2.0 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from matplotlib>=2.2->seaborn) (8.2.0)\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from matplotlib>=2.2->seaborn) (2.8.1)\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from matplotlib>=2.2->seaborn) (1.3.1)\n",
      "Requirement already satisfied: cycler>=0.10 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from matplotlib>=2.2->seaborn) (0.10.0)\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from matplotlib>=2.2->seaborn) (2.4.7)\n",
      "Requirement already satisfied: six in /home/thor/.local/lib/python3.6/site-packages (from cycler>=0.10->matplotlib>=2.2->seaborn) (1.15.0)\n",
      "Requirement already satisfied: pytz>=2017.2 in /home/thor/anaconda3/envs/iglenv/lib/python3.6/site-packages (from pandas>=0.23->seaborn) (2021.1)\n"
     ]
    }
   ],
   "source": [
    "!pip install seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "a4f7fc79",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "03d000a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD4CAYAAABPLjVeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deZxeRZX3v0cJa4CwprM3IIzIrmHNiA1hX4wsAoKoiEZBBVGMzvC+DjrDCBFhUF4EZHEkENZAEiAhQAgg+xq2KGt30p102JdIUJKc9497e6Zt61bdJ0/30083v28+95On77l1q566p05X1z2njrk7QgghasPHersBQgjxUUJGVwghaoiMrhBC1BAZXSGEqCEyukIIUUNWqUEdco8QQpTFqr3Bh6+/XNrmDNhw06rrq5RaGF0+fP3l4PkBG24KwF/n3RWUr7blHnxw7xVB2eqfPRaApVMnBuVrjJuQyS89NSw//mzeP/sbQdmap14CwF9O+2JQvtYZ1wHw3on7B+VrXzCDt4/ZMygbdOVsAN448HNB+Qa33M3ipqagDGDwnDks2HFsUDbikTsBeGW7vYPyTebezvNb7ld47y3mzWTe5gcEZVu+cCsAz2x6UFC+9cs381TjwYX33rZ5eqF82+bpAFF5Ub0ddafaXfS9t5g3M9pfQLK/i57X4Dlzos8ZSOpJTMdS+pnS79jYSI2r2LgsGs+QjemUPejv1MToCiFEzVixvLdbECVpdM3sk8A4YBjZUsFCYJq7z+vhtgkhROUsX9bbLYgSfZFmZj8GriZbZ3kYeCT/PNnMftLzzRNCiMpwX1H66A1SM93jga3c/cPOJ83sHOBZ4MxQITMbD4wHuOiiizju0L26oalCCFGCFb1jTMuSMrorgKFAS5fzQ3JZEHe/GLi448eihXMhhOh2emkGW5aU0f0+cKeZvQAsyM+NBD4BfLcnGyaEECtFX36R5u4zzWwLYCeyF2kGtAKPuHt9fzMhxEeTOp/pWg22dlRwhBCiLFUHK/z1pQdL25zVNtulfwZHCCFEzejjL9K6hQ8euSF4fvUdDwNg6c3nBOVrHPQDll4+ISw7LouYef+XXw/K1/zRZQAs+dEhQfnAX97Iu8eHI5HWuTSLRHpzXDiaaP2pWTRR++5NQXnDPXNo3j5878Yns3vHIqieGDkuKAPYYf5U7h9yWFC226Ksn+cMDkcqNS2+jlmDjyq89z6Lr+aWwV8Kyg5cPBmAqQ1HB+Xj2q9iSoEM4NCI/ND2qwCi8qJ6O+pOtbvoe++z+OpofwHJ/i56XjvMn5qMlEvpSUzHUvqZ0u/Y2EiNq9i4LBrPkI3plD2omjpfXtBMVwjRv+jLL9KEEKLPoZmuEELUkDoPA5bRFUL0L+r8RZo2MRdC9Cvcl5c+YpjZCDO7y8zmmdmzZnZy4JomM3vHzJ7Mj5+m2qeZrhCif9F9a7rLgB+6++NmtjbwmJnd7u7PdbnuXncv3vC5CzK6Qoj+RTctL7j7ImBR/vk9M5tHFpnb1ehWhCLShBD1RNURYh88dlNpm7PG6EO+Rb4jYs7F+YZdf98os0bgHmBrd3+30/km4Aay7REWAqe6+7OxOmsy01167c+D59c4Ilv+eP83Jwbla37vgmRKkne/uU9Qvs7vZgHwxsEF6VKm383C3fYIyoben6UbefFT+wbln3juNgCeHPX5oHz7lmnc13B4UDam/Xog7qyfCjK4cuiXg7JjFk4C4NLhYfnxrZO4YERYBnDigkmcNzIsP3l+du9fFch/OH8SZxfIAE6dP4mJo8LyCS3ZvYvKnzp/UmG9HXWn2l30vU9cMCnaX0Cyv2NBHbHnDCT1JKZjKf1M6XdsbKTGVWxcFo1nyMZ0yh5UzfIP09fkdNkRMYiZDSQzrN/vbHBzHgdGufsSMzsAuAnYPHY/vUgTQvQvVqwofyQwswFkBvdKd5/SVe7u77r7kvzzrcAAM9swds+k0TWzT5rZ2Nzadz5fnOFQCCF6C19R/ohgZgZcCsxz92Bss5k15NdhZjuR2dQ3YveNLi+Y2UnAd4B5wKVmdrK7T83F/wnMjLZaCCFqTff56Y4BjgWeNrMn83P/SranOO5+IXA4cIKZLQOWAkd54kVZak33m8Bn8vWKRuB6M2t09/OILHh3Tddz7KBELUII0V10n/fCH0m82HP384HzK7lvyuh+vNN6RXP+pu56MxsVa0zXdD1FC+dCCNHdeAUv0nqD1Jpuu5lt3/FDboAPAjYEtunJhgkhxErRTWu6PUVqpvsVsqiM/8HdlwFfMbOLeqxVQgixstT53gsKjhBC1BNVB0csvePC8sERe327f6brSWZ3+EHYAXzgOdN4+0thB+9BkzMH78VNTUH54DlzAHh5m7CT96ZPz4o6ngPc0xB2AN+9PQvMmN4QzlZwcPvkpEP9RQUO+d9qTQcCnDHqmKDstJYrAfg/jWFn/f9ovooJjeE2A0xsnswpjWFn/nObM2f+kxqPDMp/3XwN3y2QAZwfkZ/ffA1AVF5Ub0fdqXYXfe+JzZOj/QUk+zsWMBJ7zpAOvIjpWEo/U/odGxupcRUbl0XjGbIxnbIHVVPnM13tvSCE6F9oE3MhhKghy7SJuRBC1A7NdIUQooZoTVcIIWqIZrpCCFFDNNMVQogaUuczXQVHCCHqieqDI679efngiCN+2j+DI4QQomb0/ESyKmpidN87KZwoc+1f3wzAW4c1BeXr3TCHtl33DMqGPTAbgD9tcUBQ/snnbwXgwaGHBuW7LJzCzIJUKvvlqVQmDw1HIn1pYRaJdGFBCphvL0inpolFQf2gILoK4Jzmqzmx8Yig7ILmawH4emM4Bcxlzdfz5VHh/gCY1DKFI0d9ISi7puUmAA4viHK6vmUahxbIAKa0TOOQkQcHZTfOnw5QWH5Ky7TCejvqTrW76HtPapkS7S8g2d9Fz+uc5quT0W4pPYnpWEo/U/odGxupcRUbl0XjGbIxnbIHVaM1XSGEqCF1bnQrzpFmZn/oiYYIIUS30Je3djSzaV1PAXuY2SAAdy/+m08IIXqD5ct7uwVRUssLw4HngEvIvBAMGA38Klaoa7qe4n2thBCim6nz5YWU0R0NnAycBvzI3Z80s6XufnesUNd0Pe+d1HXCLIQQPURfNrruvgI418yuy/9fnCojhBC9Sp0HR5QyoO7eCnzRzA4E3u3ZJgkhxMrjK+rbT1cRaUKIeqLqCLH3Lzy5tM1Z89vn9c+ItDcP+Vzw/Po3ZkvDrTuHHa2HPzSb5zY7MCj71Eu3AHBfQ9ixfUx75tg+tSHsnD6u/SouHxZ2PD+uLXNMP7cgDcsp8zP56QVpXE5vuTLqMA/wzcZwqpXfNV+XDGCIBREAHDQy3Gc3z7+FfUbsV3jvWQtmssfwvYOyu1pvB2D3YWOD8nva7mTMsPBzBLivbXah/L62LNAlJi+qt6PuVLuLvvesBTOj/QXxoA2IB17EnjPEAysgrmMp/Uzpd2xspMZVbFwWjWfIxnTKHlRNH/deEEKIvkVffpEmhBB9DhldIYSoIdrwRgghaohmukIIUUPq3GVMRlcI0b+Q94IQQtQOr/PlBQVHCCHqiaqDFf5yxldK25y1TvtD/wyOEEKImtEf9l6oloW77RE8P/T+uwCYt3k4NciWL9xaVcQZwKXDw1E5x7emU+qcVpBq5Yw81cp3G48Mys9vvoZjC6KUrmiZAsSjnPYfsX9QBjBjwQyahu8VlM1pvQOAXYeF+/uBtrsYPeSzhfd+dNG9bNewW1A2t/1+ALYevEtQ/sziB9ly450K7z3v1YcL5fNefRggKi+qt6PuVLuLvveji+6N9heQ7O+i5zVjwYxkNFtKT2I6ltLPlH7HxkZqXMXGZdF4hmxMp+xB1ehFmhBC1JBl9f0iLZqux8x2NrN18s9rmNnPzGy6mZ1lZuvWpolCCFEBdZ6uJ5Uj7TLg/fzzecC6wFn5uct7sF1CCLFyrPDyRy+QWl74mLsvyz+PdvdP55//aGZPFhXqmq4nnHBZCCG6n+5yGTOzEcAfgAZgBXCxu5/X5Rojm5AeQDYZ/Zq7Px67b2qm+4yZHZd/nmtmo/OKtgA+LCrk7he7+2h3Hz1+/PhEFUII0Y1030x3GfBDd98S2AX4jpl9qss1+wOb58d44Lepm6aM7jeAz5nZS8CngAfM7GXgd7lMCCHqi24yuu6+qGPW6u7vAfOAYV0uGwf8wTMeBAaZ2ZDYfVM50t4BvmZmawOb5te3uvviaGuFEKK3qCAMuPNSaM7FeWLdrtc1AjsAD3URDQMWdPq5NT+3qLBORaQJIeqIqiPE3vv+waVtztr/NT1Zn5kNBO4GznD3KV1ktwC/cPc/5j/fCUxw98eK7lcTP92Xtt43eH6zZ24D4KGhYQfxnRdO4ZbBXwrKDlw8GYDfF6Qk+VqekuTsgpQmp86flHQuP7HxiKD8guZrgXialnEjw68Pp86/GYCxw/cJyu9snVXorA+Zw/4ODWOCsifa7wNgq8E7B+XPLn6IzTf6TOG9X3jtMTbZYLug7JU35gIwcv1tgvL5bz7NsPW2Krx321vPFsrb3noWICovqrej7lS7i773C689Fu0vINnfseCK2HMGknoS07GUfqb0OzY2UuMqNi6LxjNkYzplD6qmG70SzGwAcANwZVeDm9MKjOj083BgYeyeqTVdIYToW6xYUf6IkHsmXArMc/dzCi6bBnzFMnYB3nH3wqUFUESaEKK/0X0z3THAscDTnVxk/xUYCeDuFwK3krmLvUjmMnZc4D5/h4yuEKJ/0U1GN1+nja75evZS7DuV3FdGVwjRr/Dl2mVMCCFqh3YZE0KI2uEyukIIUUPq3OgqOEIIUU9UHRzxzrFjS9ucda+4U+l6hBCiGnyZXqTx9CYHB89v88p0AGYPDkfW7Ln4Wq4ZckxQduSiKwH4zYhw5Mz3FmSRM6ePCpc/veVKTipIhfLr5muAdCqVWDRRKhJpl6FNQfmDC+ewbcOuQRnAU+0PsMVGo4Oy5197FIBRG2wblLe88RRDBnXdJOl/WfT2c2y07j8FZa+982cA1l9786D8zfdeYL2Bnyi891tLXmTdgZsFZe8seQmgsPxbS14srLej7lS7i773orefi/YXkOzvouf1VPsD0ecM8chEiOtYSj9T+h0bG6lxFRuXReMZsjGdsgdVU982VzNdIUT/ok+/SDOzVYGjgIXufoeZHQ3sRrbF2cXuXrinrhBC9Ap9fKZ7eX7Nmmb2VWAgMAUYC+wEfLVnmyeEEJXRp2e6wDbuvq2ZrQK0AUPdfbmZTQLmFhXqmq6neIVSCCG6mT4+0/1YvsSwFrAmWWLKN4HVgAFFhfJNgDs2Avanf9FNC+RCCJHgf7I61ikpo3sp8Cfg48BpwHV5up5dgKt7uG1CCFExvZRZvTSpdD3nmtk1+eeFZvYHYC/gd+7+cC0aKIQQFVHnRlcRaUKIeqLqCLHX9v5caZuz0e1398+ItEeGHRI8v2PbjUA89cflBWlDjsvThkwcFZZPaMnkExrD957YPJlvNB4elF3SfD0AXxw1Lii/rmUqAPuO2D8ov23BDMYM2zMou69tNgDbNewWlM9tv7/QGR8yh/xYyhyAhkFbBuXtb89LBhmss9amQdm7f3kZgLXWbAzK//J+M6uvPrLw3h98ML9Q/sEH8wGi8qJ6O+pOtTsW1BHrL4inKIJ48ETsOQNJPYnpWEo/U/odGxupcRUbl0XjGbIxnbIH1dKnlxeEEKKv4ctrPnmtCBldIUS/QjNdIYSoIb5CM10hhKgZmukKIUQNcddMVwghaoZmukIIUUNW1Ln3goIjhBD1RNUWs+XTe5W2OaMev6N/BkcIIUStkPcCcP+Qw4Lnd1t0AwA3NRwdlH+h/SouGh6OfPlWaxYZ8+8FKUf+b0uWzueUxqOC8nObr+arjeF2/Xdz1q5YqhSAPYbvHZTf1Xo7Ow7dPSh7ZOE9AGw1eOeg/NnFD7HphjsEZQAvv/5ENPUMxKOvBq65SeG9l7z/SjJqbNXVhgflf/trKwNWHVZ47w//1lYo//BvbQBReVG9HXWn2l30vZe8/0q0vyCe6gcofF4vv/5E9DkDST2J6VhKP1P6HRsbqXEVG5dF4xmyMZ2yB9XS83+8V4dmukKIfoVmukIIUUPkMiaEEDVkeZ17L3wsJjSzdc3sTDP7k5m9kR/z8nODIuXGm9mjZvboxRdfXHSZEEJ0O+5W+ugNokYXuBZ4C2hy9w3cfQNgj/zcdUWF3P1idx/t7qPHjx/ffa0VQogEvsJKH71Byug2uvtZ7t7eccLd2939LKB481QhhOgl3MsfvUHK6LaY2QQzG9xxwswGm9mPgQU92zQhhKicep/pRiPSzGw94CfAOGDj/PRiYBpwpru/VaKOOveaE0LUEVVbwqc3Obi0zdnmlen1FZGWG9Uf58ffYWbHAZeXqeS+hnDakDHtWdqQKQXO1Ie2X8WFI8JO2N9ekAVH/LzAifunuRP3yQUO4Oc1X82xow4Nyq5omQLAwQXO59Nz5/Om4XsF5XNa72D0kM8GZY8uuheALTfeKSif9+rDjNpg26AMoOWNp5LpZdYb+Img/K0lLybT3qy2+oig7K8fZH/YxAIYVokERyyLyJflwRExeSrwItXuWJqhWH9BPP0RUPi8Wt54KvqcgaSexHQspZ8p/Y6NjdS4io3LovEM2ZhO2YNqqffgiNTyQoyfdVsrhBCim1jhVvroDaIzXTN7qkgEDC6QCSFEr9HXgyMGA/uSuYh1xoD7e6RFQghRBfW+vJAyujcDA939ya4CM5vTIy0SQogq6M5lAzO7DDgIeNXdtw7Im4CpwCv5qSnu/vPYPVMv0o6PyIpXy4UQopdYvqKaV1X/wO+B84E/RK65193DbzUDdGvrhBCit/EKjuS93O8B3uzO9snoCiH6FZV4L3TeJyY/Vmbfgl3NbK6ZzTCzrVIXK12PEKKeqHpB9r6Gw0vbnDHt1yfrM7NG4OaCNd11gBXuvsTMDgDOc/fwrvg5mukKIfoVKyo4qsXd33X3JfnnW4EBZrZhrIzS9QRQup5/ROl6/hGl6/lH6iJdT/WT5dKYWQOw2N3dzHYim8i+ESujTcyFEP2KZd3rMjYZaAI2NLNW4N+AAQDufiFwOHCCmS0DlgJHeWLNVkZXCNGv6M6Zrrt/KSE/n8ylrDQyukKIfkV3rNX2JKl0PeuY2S/M7AozO7qL7IJIOaXrEUL0Co6VPnqDlPfC5WQuHDcAR5nZDWa2Wi7bpaiQ0vUIIXqLWnovrAyp5YXN3L3jVeNNZnYaMNvMPt/D7RJCiJVieS/NYMuSMrqrmdnH3H0FgLufkb/BuwcY2OOtE0KICumlLDylSaXrmQjMcvc7upzfD/hNKvIiRxFpQoiyVG0ypzYcXdrmjGu/qu7S9UwoOD/TzP6zbCWPDDskeH7HthsBuGVw2CvjwMWTuXxY2An7uLYsOGLiqLB8Qksmn9AYvvfE5sl8ozGcNuSS5ixtyBdHjQvKr2uZCsC+I/YPym9bMIMxw/YMyu5rmw3Adg27BeVz2+9ni41GB2UAz7/2KCPX3yYom//m00A8vUxRIABkwQDrrLVpUPbuX14G4mlvigIUIAtSSAUwxOSpNEOpdscCIFLpeFL9XfS8nn/t0ehzBpJ6EtOxlH6m9Ds2NlLjKjYui8YzZGM6ZQ+qpd5neUrXI4ToV/TpF2lK1yOE6GussPpe1FW6HiFEv2J5bzcggdL1CCH6FfXuvaB0PUKIfsWKPu6nK4QQfYp6916Q0RVC9CvqfXlB6XqEEPVE1Sbz98O+XNrmfK1tUn0FR3QXT29ycPD8Nq9MB2D24COC8j0XX8s1Q8I72B+5KNvB/jcjwk7a31uQOXGfXrAD/uktV3JS45FB2a+brwHg2FGHBuVXtEwB4pklxg7fJyi7s3UWALsMbQrKH1w4h20bdg3KAJ5qfyDqjA8waoNtg/KWN54qzIIAWSaEjdb9p6DstXf+DMSDDNYb+InCe7+15EXWHbhZUPbOkpcACsu/teTFZFBHqt2x7A+x/oJ48ANQ+Lyean8g+pyBpJ7EdCylnyn9jo2N1LiKjcui8QzZmE7Zg2pZXuczXS0vCCH6FfW+n66MrhCiXyGjK4QQNaQbU6T1CBXvvWBmG5e4RpkjhBC9Ql/fe2H9rqeAh81sBzLPhzdD5dz9YqDD2vrTv+ieBXIhhEjR18OAXwdaupwbBjxO5goW3k9PCCF6iXr3000Z3QnAXsCP3P1pADN7xd036fGWCSHESlDvL9KSwRFmNhw4F1gA/Bsw190rmeEqOEIIUZaq56m/Glk+OOKH82sfHJF8kebure7+ReAu4HZgzR5vlRBCrCRewdEblHYZc/fpZnYHsBmAmR3n7peXKfvS1vsGz2/2zG0APDQ0HFmz88Ip0VQ+AL8vSBvytTydz9kjw/JT50/itMbwRmlnNF8FwImN4ciaC5qvBeDLBRFBk1qmRCOJIB6JtOuwPYIygAfa7mKHhjFB2RPt9wGw1eCdg/JnFz/E5ht9pvDeL7z2GJtssF1Q9sobc4F46pph621VeO+2t54tlLe99SxAVF5Ub0fdqXYXfe8XXnss2l9Asr+LntcDbXdVFXEGcR1L6WdKv2NjIzWuYuOyaDxDNqZT9qBa6n1NtyKXMXdf6u7P5D8qXY8Qou5YXsHRGyhdjxCiX7Gizl8jKV2PEKJfUe/eC0rXI4ToV9T3PFfpeoQQ/Yy+PtMVQog+xTKr77mujK4Qol9R3yZX6XqEEPVF1V62/9J4dGmb84vmq/pnup6Fu4Wdx4fefxcA8zY/ICjf8oVbua/h8KBsTPv1AExtCC8tj2vPHMAvHR528j6+dRITR4VlE1oyB/CUc/l3C9KhnN98TTKVyqGjPh+UT2mZxv4j9g/KAGYsmEHT8L2CsjmtdwBxZ/3RQz5beO9HF93Ldg27BWVz2zNnla0H7xKUP7P4QbbceKfCe8979eFC+bxXHwaIyovq7ag71e6i7/3oonuj/QUk+7voec1YMCP6nCGdEiqmYyn9TOl3bGykxlVsXBaNZ8jGdMoeVEtfdxkTQog+RX2bXBldIUQ/o969FyrOHCGEEPXMcrz0kcLMLjOzV83smQK5mdmvzexFM3vKzD6duufKpOvZoMQ1StcjhOgVujldz++B/SLy/YHN82M88NvUDaNG18zONLMN88+jzexl4CEzazGzzxWVc/eL3X20u48eP358qg1CCNFteAX/kvdyvwcIpiXLGQf8wTMeBAaZ2ZDYPVMz3QPd/fX88y+BI939E8DewK+SLRZCiBpTyUy381/l+VHpLHEYWYKHDlrzc4WkXqQNMLNV3H0ZsIa7PwLg7s+b2WoVNk4IIXqcSlzGuiTRXRlCfr7RBkSDI8zse8DBwJnA7sAgYAowFtjU3Y8t0ah69+AQQtQPVQcrnNB4RGmb89vma5P1mVkjcLO7bx2QXQTMcffJ+c9/BprcfVHR/VIb3vzGzJ4GTgC2yK/fArgJ+PdUY4UQotYsq+08bxrwXTO7GtgZeCdmcKGEn667zwHmdD1vZscBpdL1vHlI+J3b+jfeDUDrznsG5cMfms1zmx0YlH3qpVuAeGQMxCPWLi9ISXJcnpLk3IJ0JqfMz+SnjzomKD+95Up+0HhUUHZO89UAfLPxi0H575qvK0zRAlmallSU00Ejw3128/xb2GdE8YvYWQtmssfwvYOyu1pvB2D3YWOD8nva7mTMsPBzBLivbXah/L622QBReVG9HXWn2l30vWctmBntL4hHD0I8pU7sOQNJPYnpWEo/U/odGxupcRUbl0XjGbIxnbIH1VLmBVlZzGwy0ARsaGatZMl5BwC4+4XArcABwIvA+8BxqXtWExzxM0oaXSGEqBXdGRzh7uFkcP8rd+A7ldxT6XqEEP2K7pzp9gRK1yOE6FfUexiw0vUIIfoVy3t+u9qqULoeIUS/Qls7CiFEDenra7pCCNGnqPc1XaXrEULUE1VHpH1x1LjSNue6lqn9M13PeycdFDy/9q9vBuCtw5qC8vVumEPbrmFH62EPZA71f9oinBrkk8/fCsCDQ8OO67ssnMLMwWHH9P0WZ47pk4eGHdO/tPBKAC4cEXY+//aCdCqg/1OQauU/mq8qdJiHzGn+xMYjgrILmq8F4OuNYcf2y5qvTwZeHDnqC0HZNS03AXB4QaDA9S3TCoMIIAskOGTkwUHZjfOnA/EghKJ6O+pOtTsWwBDrLyDZ37EAh9hzhnRKnZiOpfQzpd+xsZEaV7FxWTSeIRvTKXtQLVpeEEKIGtKnvReEEKKvIe8FIYSoIfX+Iq1HcqQpXY8QorfozswRPUEqXc9oM7vLzCaZ2Qgzu93M3jGzR8xsh6JyStcjhOgtVuClj94gtbxwAdlWZoPI9lo4xd33NrOxuWzXHm6fEEJURA3cYKsitbwwwN1n5Luiu7tfT/bhTmD1Hm+dEEJUSHemYO8JUul6HiCb6a4LnA2c7O435ZmAf+Xuo0vUUd+/doQQ9UTVwQp7jdi3tM25Y8FtdRcc8W1gItkLwX2BE8zs90Ab8M2ebZoQQlROvS8vpHYZm0tmbDs4OT860vWU2lP3/V9+PXh+zR9dBsCSH4SjjQaeM423v7RHUDZo8l0ALG5qCsoHz5kDwMvb7BOUb/r0LJ4siHLaPk/Dck9DONXK7u1ZqpXpDeFN5Q9un8yVQ8ORRMcszCKNLhoeln+rdRK/KkjDAvDD+ZM4oyCFy2ktWSRSLApqQmPxRvgTmydzSkF01bl5+piTGo8Myn/dfA3fLZABnB+Rn998DUBUXlRvR92pdhd974nNk5NRY6n+LnpeP5w/KfqcgaSexHQspZ8p/Y6NjdS4io3LovEM2ZhO2YNqqXc/3Wpcxn7Wba0QQohuot5dxpSuRwjRr+jrYcBK1yOE6FPU+/KC0vUIIfoVfdroKl2PEKKv0ae9F4QQoq/Rp2e6QgjR16j3TcyVrkcIUU9UHSH26SH/XNrmPL7oj3UXkdYtLL3258HzaxzxUwDe/82JQfma37uAv5wWdgBf64zMAfzdb4YdvNf53SwA3jj4c0H5BtPvZuFuYQfvofdngRcvfmrfoB4CpDoAAAvSSURBVPwTz90GxJ3P72sIp4AZ056lgJlVkEpln8VXM6WheLn80Parkg71lxY45B/fOokLCtK/AJy4YBLnFTj6nzw/u3csEODsSFDHqfPTKYyKyp86Px0wkmp30fc+ccGkaH9BOoCh6Hkd2n5V9DkDST2J6VhKP1P6HRsbqXEVG5dF4xmyMZ2yB9WiNV0hhKghWtMVQogaUu9rujK6Qoh+xYo6X15IZY5Y18zONLM/mdkb+TEvPzcoUk7peoQQvUKf3nsBuBaYDTS5ezuAmTUAXwWuA/YOFXL3i4EOa+tFC+dCCNHdLPf6Tk2Z2mWs0d3P6jC4AO7e7u5nASN7tmlCCFE5K9xLH71Byui2mNkEM/ufHcXMbLCZ/RhY0LNNE0KIyqn35YVUup71gJ8A48h2HHNgMTANOMvd3yxRR32vagsh6omqgxU22/DTpW3OS68/Xl/BEe7+lpldDtwOPOjuSzpkZrYfMLOH2yeEEBXRp13GzOwk4DvAPOASMzvZ3afm4v+kpNH94JEbgudX3/EwAJbefE5QvsZBP2Dp5RPCsuMmAiVSAf3okKB84C9v5N3jg+8BWefS2wF4c1w4Ymf9qXcD0L57U1DecM8cmrcP37vxyeze8zY/ICjf8oVbeWLkuKAMYIf5U7l/yGFB2W6Lsn6eMzgcLdS0+LrCCCnIoqRuGRxOD3Pg4skATC2IvhrXflUyki4WuQXxyK6iejvqTrU7FhkW6y8g2d9Fz2uH+VOjzxlI6klMx1L6mdLv2NhIjavYuCwaz5CN6ZQ9qJblvrxb7tNTpLwXvgl8xt2XmFkjcL2ZNbr7eXTDnwFCCNHd1HsYcOpF2sc7lhTcvRloAvY3s3OQ0RVC1CEr8NJHCjPbz8z+bGYvmtlPAvImM3vHzJ7Mj+QGEqmZbruZbd+ROSKf8R4EXAZsk2yxEELUmO6a6ZrZx4H/RxaP0Ao8YmbT3P25Lpfe6+4Hlb1vaqb7FaC98wl3X+buXwF2L1uJEELUim70090JeNHdX3b3vwFXk3lyVUXU6Lp7a+fAiC6y+6qtXAghuptK/HQ7b1mQH+M73WoYfx+P0Jqf68quZjbXzGaY2Vap9mnDGyFEv6KSMOAuWxZ0JfTequv0+HFgVL70egBwE7B5rM7U8oIQQvQp3L30kaAVGNHp5+HAwi51vdvJ2eBWYICZbRi7qdL1CCHqiaq9otZfe/PSNufN914orM/MVgGeB8YCbcAjwNHu/mynaxqAxe7uZrYTcD3ZzLewDTVZXvjw9ZeD5wdsuCkAf513V1C+2pZ78MG9VwRlq3/2WACWTp0YlK8xLnPeXnrpqWH58Wfz/tnfCMrWPPUSIJ6SBOC9E/cPyte+YAZvH7NnUDboytkAvHFgQaqUW+5mcVNTUAYweM4cFuw4Nigb8cidALyyXdgpfpO5t/P8lvsV3nuLeTOTzvzPbBp+Sbv1yzfzVOPBhffetnl6oXzb5ukAUXlRvR11p9pd9L23mDcz2l9Asr+LntfgOXOizxlI6klMx1L6mdLv2NhIjavYuCwaz5CN6ZQ9qJbumki6+zIz+y5wG/Bx4DJ3f9bMvp3LLwQOB04ws2XAUuComMEFrekKIfoZ3ZmuJ18yuLXLuQs7fT4fOL+Se8roCiH6FfUekSajK4ToV/TpTczNbB0z+4WZXWFmR3eRXRApp3Q9Qoheod43MU/NdC8HXgBuAL5uZoeRvb37K7BLUaGu6XqKFs6FEKK76evLC5u5e8d+azeZ2WnAbDP7fA+3SwghVoo+vZ8usJqZfcw9WyRx9zPMrBW4BxjY460TQogKqfeZbipdz0Rglrvf0eX8fsBv3D0a7pZT3z0ghKgnqg6OWGXVYaVtzrK/tdV8i9rUhjcTgFYzG2tmAzudnwmcVLIO63yY2be6nit7VFO2N+v+KLb7o/id1e5uKVs1y/7WZmWP7qivYhJxyd8D/ky2iUMzMK6T7PFKYpw7lXt0ZcpVW7Y36/4otvuj+J3V7tq3uy8eqTXd8ShdjxBCdBspo/t36XrMrInM8I5CRlcIISomtbVju5lt3/FDboAPAjZk5dP1VBMtUW2kRW/V/VFs90fxO/dm3R/Vdvc5Ut4Lw4FlHsgeYWZjXNkjhBCiImqxn64QQogcZY4QQogaUjOjm8ofnyg7wszuMrN5ZvasmZ28EvV/3MyeMLObV6LsIDO73sz+lLdh1wrKnpK3+Rkzm2xmq0euvczMXjWzZzqdW9/MbjezF/L/16uw/C/zdj9lZjea2aCyZTvJTjUzL0pDUlTWzL6XP/Nn80CbStq9vZk9aGZP5psn7RQoF9SLsn0WKZ/ss5ROxvosVrZMn0XaXabPVjezhy1LpPismf2sbJ9FypbVsWD5Mn3Wr6iFXxrZrusvAZsCqwJzgU9VUH4I8On889pkKTRKl8/L/QC4Crh5Jdr/38A38s+rAoNKlhsGvAKskf98LfC1yPW7A58Gnul0biLwk/zzT4CzKiy/D7BK/vmsovKhsvn5EWQ757cAG1ZQ7x7AHcBq+c8bV9juWcD++ecDgDll9aJsn0XKJ/ssppOpPovUW6rPIuXL9JkBA/PPA4CHyDavSvZZpGxZHQuWL6tn/eWo1Uy3qvzx7r7I3R/PP78HzCOcCjmIZS8EDwQuqajVWdl1yIzCpXn9f3P3tyu4xSrAGpblW1qTLontOuPu9wBvdjk9jszok///hUrKu/ssd1+W//ggWXK9snUDnAtMIBLOXVD2BOBMz3akw91frbC8A+vkn9cl0G8RvSjVZ0Xly/RZQiejfRYpW6rPIuXL9Jl77gZKZvgG5OWSfVZUtgIdK6obSuhZf6FWRrds/vgklgVp7ED2W7Is/0X2QFdmd+NNgdeAy/PliUvMbK0yBd29DTgbmA8sAt5x91kV1j/Y3Rfl91sEbFxh+c58HZhR9mLLdpNrc/e5K1HXFsBnzewhM7vbzHassPz3gV+a2QKyPvyXRFsb+V+9qLjPInqV7LPOZSvtsy71VtxnXcqX6jPLltqeBF4Fbnf30n1WULYz0f4Kla9Sz/octTK6oUCKin+jWbb/ww3A99393ZJlDgJedffHKq0vZxWyP31/6+47AH8h+/OrTN3rkc0gNgGGAmuZ2ZdXsh1VYdm2nMuAK0tevyZwGvDTlaxyFWA9sj8/fwRca2aVBNScAJzi7iOAU8j/0ihoa8V6UaZ8mT7rXDa/tnSfBeqtqM8C5Uv1mbsvd/ftyWakO5nZ1mXamypbpr8C5belOj3rc9TK6Cbzx6cwswFkCnalu0+poOgY4PNm1ky2rLGnmU2qoHwr0NrpN/r1ZEa4DHsBr7j7a+7+ITAF2K2CugEWm9kQgPz/wj/TizCzr5IFtRzj7mV/2W1G9stibt53w4HHLUs5XYZWYEr+J+XDZH9lVPKC5Ktk/QVwHdkS1T9QoBel+6xIr8r0WaBs6T4rqLd0nxWUL9VnHeTLZHOA/ahQz7qUrVjHOpXvmJSsrJ71OWpldB8BNjezTcxsVeAoYFrZwvlv+0uBee5+TiUVu/u/uPtwd2/M653t7qVnm54Fhiwws3/KT40FnitZfD6wi5mtmX+HsWTrb5UwjWwwkf8/tZLClm3D+WPg8+7+ftly7v60u2/s7o1537WSvbz5h0CZAm4C9szbsAXZC8jXK2j6QqAjf/meZBlM/o6IXpTqs6LyZfosVLZsn0XaXarPIuXL9NlGHd4FZrYG2cTgT5Tos6KyZXWsoPwTVepZ38Nr9MaO7G3q82ReDKdVWPafyZYjngKezI8DVqINTayc98L2wKN5/TcB61VQ9mdkSv0McAX5m+mCayeTrf1+SKZ8xwMbAHeSDaA7gfUrLP8i2Xp6R79dWLZsF3kzxd4LoXpXBSbl3/txYM8K2/3PwGNkni4PkW28VEovyvZZpHyyz8roZFGfReot1WeR8mX6bFvgibzsM8BP8/PJPouULatjwfJl9ay/HIpIE0KIGqKINCGEqCEyukIIUUNkdIUQoobI6AohRA2R0RVCiBoioyuEEDVERlcIIWrI/wf8FAvCMCiI5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = sns.heatmap(ellip_base[0], linewidth=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bac2c608",
   "metadata": {},
   "source": [
    "- break bounding box into grid\n",
    "- record every voxel center\n",
    "- compute distance voxel center from mesh\n",
    "- marching on voxel grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "11defe09",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.10423391,  0.03628614,  0.06471388],\n",
       "       [ 0.122075  ,  0.03628614,  0.06471388],\n",
       "       [ 0.122075  ,  0.03628614, -0.06815754],\n",
       "       [-0.10423391,  0.03628614, -0.06815754],\n",
       "       [-0.10423391,  0.2060993 ,  0.06471388],\n",
       "       [ 0.122075  ,  0.2060993 ,  0.06471388],\n",
       "       [ 0.122075  ,  0.2060993 , -0.06815754],\n",
       "       [-0.10423391,  0.2060993 , -0.06815754]])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bbox_verts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "13f5fa05",
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3ad0f6b",
   "metadata": {},
   "source": [
    "### Marching Cube method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "44885adb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import igl\n",
    "from meshplot import plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "d9998a53",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cube(x,y,z,delta):\n",
    "#     print(x,y,z,delta)\n",
    "    return np.array([[x+delta, y+delta, z], [x+delta, y, z], [x, y, z], [x, y+delta,z],\\\n",
    "        [x, y+delta, z+delta], [x,y,z+delta], [x+delta, y+delta, z+delta],[x+delta,y,z+delta]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "87fbadd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# grid coords (x+delx, y+dely, z), (back top right)\n",
    "# (x+delx, y, z) back bottom right\n",
    "# (x, y, z) back bottom left\n",
    "# (x, y+dely, z) back left top\n",
    "\n",
    "# (x, y+dely, z+delz) front top left\n",
    "# (x, y, z+delz) front bottom left\n",
    "# (x+delx, y, z+delz) front bottom right\n",
    "\n",
    "def create_grid( bounds, delta=.1):\n",
    "    min_bounds = bounds[0]\n",
    "    max_bounds = bounds[1]\n",
    "    x_b = round((max_bounds[0] - min_bounds[0])/delta+1)\n",
    "    y_b = round((max_bounds[1] - min_bounds[1])/delta+1)\n",
    "    z_b = round((max_bounds[2] - min_bounds[2])/delta+1)\n",
    "    grid_cells = []\n",
    "    centers = []\n",
    "    points  = []\n",
    "\n",
    "#     print(bounds)\n",
    "    for i in range(x_b):\n",
    "        for j in range(y_b):\n",
    "            for k in range(z_b):\n",
    "                x = min_bounds[0] + i * delta\n",
    "                y = min_bounds[1] + j * delta\n",
    "                z = min_bounds[2] + k * delta\n",
    "                points.append([x,y,z])\n",
    "                grid_cube = get_cube(x,y,z,delta)\n",
    "                center = grid_cube.mean(0)\n",
    "                centers.append(center)\n",
    "                grid_cells.append(grid_cube)\n",
    "    \n",
    "    return np.array(grid_cells), np.array(centers), np.array(points),[x_b,y_b,z_b]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "0f68adff",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9261, 8, 3) (9261, 3) (9261, 3)\n"
     ]
    }
   ],
   "source": [
    "bounds = np.array([[0., 0., 0.], [2., 2., 2.]])\n",
    "grid_cells,centers,points,sizes = create_grid(bounds)\n",
    "print(grid_cells.shape, centers.shape, points.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "30ed32b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "v,t,f = igl.read_mesh('bunny.mesh')\n",
    "\n",
    "zmin = v[:,-1].min() + 0.1 * v[:,-1].min()\n",
    "zmax = v[:,-1].max() + 0.1 * v[:,-1].max()\n",
    "ymin = v[:,-2].min() + 0.1 * v[:,-2].min() \n",
    "ymax = v[:, -2].max() + 0.1 * v[:, -2].max() \n",
    "xmin = v[:, -3].min() + 0.1 * v[:, -3].min() \n",
    "xmax = v[:, -3].max() + v[:, -3].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fa3ed17",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "8f2c6346",
   "metadata": {},
   "outputs": [],
   "source": [
    "bounds = np.array([[xmin, ymin, zmin], [xmax, ymax, zmax]]) # bounds for bunny\n",
    "grid_cells,centers,points,sizes = create_grid(bounds,delta=0.006)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "76bbb528",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((26013, 3), (26013, 8, 3), (26013, 3))"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "centers.shape, grid_cells.shape, points.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "93283481",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = igl.signed_distance(centers, v, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "198e4927",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(26013, 3)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "centers.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "5a8a42ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "c_reshape = d[0].reshape((sizes[0],sizes[1],sizes[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "4dd4be3b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(39, 29, 23)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_reshape.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "3f97f3d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fb178a96438>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAG+CAYAAABcRDoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dXaiu6Vkf8P/1vGvtPckk6kRNGmJarYRSkTaWIRRSSopVoifRgsUcSArCeGBKBA8aPKhpoSDFj54UIZJgWvxAiNYchLYhWFKhWCchzUdHmyBRY4ZMJZR8zl7rfZ+7B3tFxjhrz76vZ+211p79+8Gw1373vua+n/u5n+f9r2etva4aYwQAgLu3XPUEAADuNwIUAMAkAQoAYJIABQAwSYACAJgkQAEATDq6zMF2L354HH3jI5c5ZHK//JSGURtqL7duw0wvfa4PjE0n5fLG27LN28dYV7DRl96Y1RxzaY6XJEutrbpdc113S2+822P2ao+adUvzxtOdZ5Ls0j0fzbrmeFvGPGqva/+i3DWfF33wI7f+fIzxzc/2Z5caoI6+8ZH8tX/5z+cL1w13sm7tZYeS0/7DwOW0Oea+tzYb7g2pQ69uuYK5Xnpo25Khd726ddc7yNG8c6wb7jjrjeZcbzY3wXF/8yw3exv9xs19q+6hG82bQJKHb5606r7u5tO9uhu9uiR5yY0vt+q+4ahX96LdrVbdi3f9Y/z63ZdadV/XHPMblt7a3K79SqvuJbvefv36pXmjS/L1ywtadbuXf/KPz/szX8IDAJi0KUBV1eur6g+r6pNV9daLmhQAwHXWDlBVtUvy75N8X5LvSPLGqvqOi5oYAMB1teUJ1GuSfHKM8UdjjJMkv57kDRczLQCA62tLgHpFkj99xu8/ffYaAMDz2pYA9Wz/Zuiv/DOZqnqsqh6vqscPX+z96wIAgOtkS4D6dJJXPuP335LkM1/7l8YYbx9jPDrGeHT3ooc3DAcAcD1sCVC/n+RVVfVtVXUjyQ8nec/FTAsA4Ppq/1i7Mca+qt6c5L8k2SV55xjj4xc2MwCAa2rTTyIfY7w3yXsvaC4AAPcFP4kcAGCSAAUAMOlSmwlnVOpWoxlgswFtklS3mXCzj2i78XuzWW6ypUFvc7wrmGu3KfCmxsfNczm6y7OlmXD3GjnqDbo2Oy0vV3GMzf06Nkx2NO87h0Ozrr3pkv3a+zz61qH39nFy6DeE/fzpQ6267vrcana/Pu12995Qe9rs8H3Y3U/PUfpNs5Ne4+M7uZ9WDgDgWhCgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYNLRpY62JstXarqsDvM1f1E7moXrJY/XrUuynPTWZ9n3xqvm2iRJHS65bsO6ts9Jc7uODZ/OrM0ruT1k9Q5yXfonpHsfGGtzzLV/3+nWrmvvjOz3u1ZdkpzuehfX6aE35hdOH2rVJcl+9MbcNy+u/dob77Dh2cRp8xgP3RvPA+P0wv+PnkABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADApGYP954aye7WfMfo2ve7TNd6uXXpN5tvW/a9umrWdcdLkuo1fr+vzkf305Jm4/ctQ2ZtXlpXcT7613LvIMfan+xo3rPGrld3OPQ/F751etyqG6M315NDf6OfNmu7Y+6Pe+t6OvrnY3+04UbAuY7bN63zeQIFADBJgAIAmCRAAQBMEqAAACYJUAAAkwQoAIBJAhQAwCQBCgBgkgAFADBJgAIAmCRAAQBMEqAAACYJUAAAk44udbQ12X1lvoN3HfpD9ju4NwdsNnCvq+hSv+/VLfv+ZLvn8h400r6LQXtla7OZ+tIdMMna/VSoeyq7+3zDebzsfV67/vkYh17tOPRO5P6kVZYkWZtjHo57daeH/g39tLs+zQtkjN55XJt195td80awa17Mx92LeWPteTyBAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADApKPLHKzWZHdrvm7Zbxjz0KxbmwN2667Ash/Nuv6Y7XXtTXWT0fz0oqo53hUcY/d8tK+rfXNxktTSW6BaemN2jzFJsu9tnjF6cx27/ufC61FvE6xrb66HQ38PdMdsr2urKlmz4Ri7N56mpX1TTnbN2u6Yy4Y32BubLuhn5wkUAMAkAQoAYNKmL+FV1aeSfCHJIcl+jPHoRUwKAOA6u4jvgfpHY4w/v4D/DwDAfcGX8AAAJm0NUCPJf62qD1bVYxcxIQCA627rl/BeO8b4TFW9NMn7quoPxhgfeOZfOAtWjyXJ8dc9snE4AICrt+kJ1BjjM2e/PpXkt5K85ln+ztvHGI+OMR7dveDhLcMBAFwL7QBVVQ9X1Yu/+nGS703ysYuaGADAdbXlS3gvS/JbdfvHLh8l+dUxxn++kFkBAFxj7QA1xvijJH/3AucCAHBf8GMMAAAmCVAAAJMu4ieR37Uaye7WfN1y2m9T3+4a363rN4u+dN25LocN5+OS16fZhD1Jsh71iqu5PLVhXZfT5oF2h+wOt+FTtjo0z8dR7yCXW/3NM/bNuuaeG81jTJLRXNe1WTeON6zr2ttAo3kj6NatW248TUvz5rrr3rCSLM3a4+YbbLcuSW4c+rXn8QQKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJh1d6mhrsnt6vnvzctofcml2RV8OzQ7VzbINDbEvXXXXZtOgza7oG3Z4s7l5f79uWNa1uYFqaXaNb65Nd02T/rVca+8Yx65/QsauOWbzGLfs83HUvLZuNMdrno8kOYzmurZHvHxL81pemkfZHS9JjpoX9HEdLrVua+15PIECAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASRt6eM+rkexuzdftTvrdopd9s0N1s3Fzrc253k/twjcYS7Pz+1FvgerQ7/xe3S7l1Rxzwx5oNkVvj9ldmuW0V5ck2ffKRvPTxO5eTZKxa9Y178h13KtLkrXbpH7tLey40d/oo1navTxOmnVbdO87u27d0l2d5KhZe9x8g725NG8CSY6ru9HP5wkUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmHV3qaCNZ9mO6bHcyX/MXtbd6tZ15Jkmtzbl2667AetzP3aO548ZSzbreeEmS3pD9MZvjbRmzmttu93SvbpNLPh9j16tLkrW5z9dDr64O/c1Tx8173eiNueVWtzY/5+8OuTbrTpt1SbIsvdl+pVm3W7pHmRw1a28s+0utS5KHli1n5dl5AgUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTmj3D+zrd36vfLDrLaa/46OleW/TaN/t+jw0tyruq10293w87ORz1WtyPZtQfveGSJOtRs9t886rqHmOStBuNn/T2Xec63mr0Tkd7D2zZO9XcO0tz7xxu9OqSJGt3YTeMecnW5rOC7jV5WPqLc3LS23jL0nyv2/VjwFFzzBtLb8PeXHrvy0nywuWkXXseT6AAACYJUAAAk54zQFXVO6vqqar62DNee0lVva+qPnH26yP3dpoAANfH3TyB+uUkr/+a196a5P1jjFclef/Z7wEAHgjPGaDGGB9I8rmvefkNSd519vG7kvzABc8LAODa6n4P1MvGGE8mydmvL724KQEAXG/3/JvIq+qxqnq8qh7fP/2lez0cAMA91w1Qn62qlyfJ2a9PnfcXxxhvH2M8OsZ49Oihh5vDAQBcH90A9Z4kbzr7+E1JfvtipgMAcP3dzY8x+LUk/yPJ36qqT1fVjyb5mSTfU1WfSPI9Z78HAHggPOfPcB9jvPGcP/ruC54LAMB9wU8iBwCYJEABAEzqt2HuuuQO3mPpdRof1eymfrrv1d3q1W1xeLjXEbuuogt783ysR81O80lGs3Q57dVVr7H5WW3vpLTP5ZXsgV7Zuuvund54SbIc9RbocNw9yF5ZklTvNrDBhmuyO2LzUUH7/WPpP5s4LLtW3elpb3VuLf3Nc9ysfXp33Kq7tZ606pLky+vFb3RPoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMOnoqidwr9Xo1vUKa23WnR5adbeLq1e2X3t1h+aiJql9r3Zp1h093SpLkozmpxfVW9b23tkyZi75+thiLL19vu56c12PeuNtqu3eP7rn/3b1pZa169K/JtPcO2PpDdieZ5LR3K+HXW/Q06Ndqy5JTg692luHXvR4+nDcqkuSr2yoPY8nUAAAkwQoAIBJAhQAwCQBCgBgkgAFADBJgAIAmCRAAQBMEqAAACYJUAAAkwQoAIBJAhQAwCQBCgBgkgAFADCp1xK5qUayHObrln2/8/ty2mtTXifNutPGASbJvlmXJNXrNF7NtemuaZIs3S71zbKl2YU9SUazSXk1l6cO/X3eHTOjN2b1p9o2mueyex4PN/qfX67Nxu/ddT2M/j5vD9q872yZancPZOkd42hugXW34b5z2hv0cNSr2++bF0iSk12v9vSoV/f0oXlhJTlZLz7ueAIFADBJgAIAmCRAAQBMEqAAACYJUAAAkwQoAIBJAhQAwCQBCgBgkgAFADBJgAIAmCRAAQBMEqAAACYJUAAAky6+PfG9sKXze7fReLNLfbe7fbsuSe0PvbrT3ulfTvq5e9fspr7su93UN3RFbx5mu7n9YcMe6NY29133GLfonsv1uHciq3dZJUkOa3Pfjf5+7RrVG7N7fXTrNo256x1j7Xrj1b5XlyTjqHk+9r3FOez6J2TfvLZuHXoLe9KsS5KvHI7btefxBAoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJOOrnoC91z1ysauly3Hjd6S1n5t1SVJDie9MU9Oe+PtmouaZDdGq24szTG7dVvG7B1i6tDfA+3aZlmt3cIN5+Ood03WfterW3t1SZLR/dy0WbdlXZtDjt0lX8tJRvMda+wvebxDry5J6tBbn9GsW9f+c5TDoVe7P/Suraf3x626JPnisuE99hyeQAEATBKgAAAmPWeAqqp3VtVTVfWxZ7z2tqr6s6r68Nl/339vpwkAcH3czROoX07y+md5/RfGGK8++++9FzstAIDr6zkD1BjjA0k+dwlzAQC4L2z5Hqg3V9VHzr7E98iFzQgA4JrrBqhfTPLtSV6d5MkkP3feX6yqx6rq8ap6/PTWl5rDAQBcH60ANcb47BjjMMZYk/xSktfc4e++fYzx6Bjj0eObD3fnCQBwbbQCVFW9/Bm//cEkHzvv7wIAPN88589ZrapfS/K6JN9UVZ9O8tNJXldVr87tn7n8qSQ/dg/nCABwrTxngBpjvPFZXn7HPZgLAMB9wU8iBwCYJEABAExq9pruGZV0Gpyvx/3u3Ydup+leo/Fk7RVWs9N8kqTZ3bxOmy3KN3R+r/3Fd8S+k7Hrz3XLcbaG27I2+2b799Hc6Gtzrlv2znHvdlU3mnXdtUnSvbWO5rU8lv5cl8Y9+XZdb67rhned6t6ymsdYzcuq9hvuHYfmuTw0z0ezLklOT3sL+3Rzv1b19/luufj3Hk+gAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYNKGvtgNlazH852f1/2GDswn7dL7Rq3N9TnttTbf1KX+tJnZm2NW9TuNj90lz/WwoVv4vtk2fm2Oeei2qd/Qpf7mjV7doVe3Ze8sS692t+vVjWZdkizNd4Hl0Btz2XA/X5tjVnObV3e8DffIWpt7YN+sO+k/RznUrlV3uvTW52jp3yNv7S4+7ngCBQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASUeXOtpIah3TZbVuGLNZuxx6hcvpoVVX+w0HuTZr9/vLrUuSqn7tJaul+flF83yMfW/v3B6zWXvaPJdL8zzudr26pL2u1dxzy27D55e73pjjuDfmsu/PdTn05lqH+Xt5ktTavwcszW3evUV233tqwy1yudVdn94e2HA6sja33X7X2zsnu/775PHRhnvPOTyBAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgElHVz0Bzoxed+otteP0tDfeumGuXcuGluFd3eMcvY7h49DvNN5uN9916Nb151mH5qBHzdvchu7tte+NWfvenqsN12S3ttbeNVlbtmrzMKu5dZbmLTLZcr9q3s+b98ja9ec6mrXrvvfs5nDUf+Zzuu9fz+fxBAoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmNduUX7Itja2btd3O1qPZLXrsNnSKPm52fr9xo1U3Tk5adUmSQ7Mt+noFWb+5B9I8xKz9NvVjv2/X3je618jN3trUYcP56J7LMXp1a7MuSZql1R1yw1R3J93i3rXcPcax4T1r6b737HqTHfv+ZMdRc66HXt166L8PHLaclHN4AgUAMEmAAgCYJEABAEx6zgBVVa+sqt+pqieq6uNV9Zaz119SVe+rqk+c/frIvZ8uAMDVu5snUPskPznG+NtJ/n6SH6+q70jy1iTvH2O8Ksn7z34PAPC895wBaozx5BjjQ2cffyHJE0lekeQNSd519tfeleQH7tUkAQCuk6nvgaqqb03yXUl+L8nLxhhPJrdDVpKXnlPzWFU9XlWP7299adtsAQCugbsOUFX1oiTvTvITY4zP323dGOPtY4xHxxiPHt18uDNHAIBr5a4CVFUd53Z4+pUxxm+evfzZqnr52Z+/PMlT92aKAADXy938K7xK8o4kT4wxfv4Zf/SeJG86+/hNSX774qcHAHD93E0PkNcm+ZEkH62qD5+99lNJfibJb1TVjyb5kyQ/dG+mCABwvTxngBpj/G7ObyT03Rc7HQCA689PIgcAmHQ3X8K7OJWsu/mOyGuzCXuSjMZ4STKaHbHXo14mrZv9g6z1uFfY7OBeRxu2zX7fKhv7Q2+8tVmXJNXs3r3rnctqno8kyaF3nOutW73xxoa5NtVDD/UKt6zrZWtOtTYcYq2XW3f0dH+y63HzmmwvbG+8Zct7VvOxRjVvdd26JKm1+f7arFtH9/wnh8PFPy/yBAoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJOOLnW0SsZuvqxT81XrcbXqDjd72bLW0apbDxtOxdorW6q3Nrl53KtLkv2hVVan+954t056dUnG2lzY5h7IrlmXJLveRVLNurFvno/unrsKY8P5aG6d6o65Yard2qMv9w6ye09OkkN3eZbemGvzmlwO/WPs3naqXbfhmuze65pjjg1zHePi7z2eQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMCko8scbCRZGyMeNnTvrkO3rjfmYe11t9+kuzxHzfzcbReepPbN2tPeVq2l/zlCfeXpZmGvQ/nYcDnWobfR64UvbNWN0ezC3pznlegeY5JqXiN16I25nPavyaMvNwub953D2r8mx65X272fV3dZ+6eje/tIrc0Tsvb3edpj9spGd7wkY/Rrz+MJFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEzqt3/vWJL1xnxH5NrQLXptNn8/jGa2bLbvHrXrjZdkt/S6TK83euva7RifJLXvrc9ysu/VjQ1zPTQ3z75XV7v+HsiN435tQ63Nfd6sS9LvGt+8PrJh7+TQO85q7vMNOyfLafP+cbM36uHGhs/bm6ekmueyRm9tasPW6R4jd7blcj6PJ1AAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJh0dJmDjSRrY8Q67o+5H9Wqa5ZlVC+T7mr0BkwyjnqTrUNvzG5dkuxuNec6emOO412rLklyo7nxbt7oj3m/aJ6P2h8ueCJ3M2jzYt5iXVtlddpcn+b5SJLxUG+fj6W5rltOR/u+fNnvA726JNvW53muNrxP3ovbgCdQAACTBCgAgEnPGaCq6pVV9TtV9URVfbyq3nL2+tuq6s+q6sNn/33/vZ8uAMDVu5vvSNon+ckxxoeq6sVJPlhV7zv7s18YY/zsvZseAMD185wBaozxZJInzz7+QlU9keQV93piAADX1dT3QFXVtyb5riS/d/bSm6vqI1X1zqp65ILnBgBwLd11gKqqFyV5d5KfGGN8PskvJvn2JK/O7SdUP3dO3WNV9XhVPX74ypcuYMoAAFfrrgJUVR3ndnj6lTHGbybJGOOzY4zDGGNN8ktJXvNstWOMt48xHh1jPLp7wcMXNW8AgCtzN/8Kr5K8I8kTY4yff8brL3/GX/vBJB+7+OkBAFw/d/Ov8F6b5EeSfLSqPnz22k8leWNVvTq3f8D4p5L82D2ZIQDANXM3/wrvd/PsP1z+vRc/HQCA689PIgcAmCRAAQBMupvvgbo4S7K2GtX32yj3u2J3uz43u35viLLLvlvXm+vuVq/TfNLv4N7tpr7FePihXt0VzPWy1WheH6eH9pjtdb2PPk0cx7tLrUuStVm7HvcWdj3uXx/rrnn/aC5P+768Yc9137Pa73VbblfVvA88T26R99GtBQDgehCgAAAmCVAAAJMEKACASQIUAMAkAQoAYJIABQAwSYACAJgkQAEATBKgAAAmCVAAAJMEKACASQIUAMCko8sesN3dumk59LpFL/vmePveeLuTZlfrDbXVXpv+XJfD2q7tODx8s1+8NDu/30+dxqs32Wqex+V016pLkrHrno/LPyE1etfIetS7QY7j/o11bdYebvTq1qP++Vib71hrc9t191x3vKT/Hjl2vT3XrUvSfwSzNN+zNuSHZbn49x5PoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMOno0kcc8yW19ofr1i6HxkSTLPte3e6kV5cky0nvIJfT3pi19ue67pqZ/YWXn/VHNet2zcJq1l2B5bS359Yb/b0zlu66dgds1m0wjnqTXZt1SbIeN8ds1h1uXMVcm+M13yHHrleXJOOSx9w0113zImnWLUs/DCzLxV/QnkABAEwSoAAAJglQAACTBCgAgEkCFADAJAEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEkCFADApBrj8lqOV9X/TfLH5/zxNyX580ubzP3H+pzP2tyZ9Tmftbkz63M+a3Nnz5f1+RtjjG9+tj+41AB1J1X1+Bjj0auex3Vlfc5nbe7M+pzP2tyZ9TmftbmzB2F9fAkPAGCSAAUAMOk6Bai3X/UErjnrcz5rc2fW53zW5s6sz/mszZ0979fn2nwPFADA/eI6PYECALgvXIsAVVWvr6o/rKpPVtVbr3o+10lVfaqqPlpVH66qx696Pletqt5ZVU9V1cee8dpLqup9VfWJs18fuco5XpVz1uZtVfVnZ/vnw1X1/Vc5x6tUVa+sqt+pqieq6uNV9Zaz1x/4/XOHtbF/klTVQ1X1P6vqf52tz786e93eOX9tnvd758q/hFdVuyT/J8n3JPl0kt9P8sYxxv++0oldE1X1qSSPjjGeDz9PY7Oq+odJvpjkP4wxvvPstX+b5HNjjJ85C+CPjDH+xVXO8yqcszZvS/LFMcbPXuXcroOqenmSl48xPlRVL07ywSQ/kOSf5QHfP3dYm38a+ydVVUkeHmN8saqOk/xukrck+Sexd85bm9fneb53rsMTqNck+eQY44/GGCdJfj3JG654TlxTY4wPJPnc17z8hiTvOvv4Xbl943/gnLM2nBljPDnG+NDZx19I8kSSV8T+udPakGTc9sWz3x6f/Tdi79xpbZ73rkOAekWSP33G7z8dF+4zjST/tao+WFWPXfVkrqmXjTGeTG6/ESR56RXP57p5c1V95OxLfA/clxieTVV9a5LvSvJ7sX/+kq9Zm8T+SXL7qyVV9eEkTyV53xjD3jlzztokz/O9cx0CVD3Law9Eer1Lrx1j/L0k35fkx8++TAN36xeTfHuSVyd5MsnPXe10rl5VvSjJu5P8xBjj81c9n+vkWdbG/jkzxjiMMV6d5FuSvKaqvvOq53RdnLM2z/u9cx0C1KeTvPIZv/+WJJ+5orlcO2OMz5z9+lSS38rtL3nyl3327Hs4vvq9HE9d8XyujTHGZ89ubmuSX8oDvn/Ovkfj3Ul+ZYzxm2cv2z959rWxf/6qMcb/S/Lfcvt7fOydZ3jm2jwIe+c6BKjfT/Kqqvq2qrqR5IeTvOeK53QtVNXDZ9/Qmap6OMn3JvnYnaseSO9J8qazj9+U5LevcC7Xyldv7md+MA/w/jn7Ztd3JHlijPHzz/ijB37/nLc29s9tVfXNVfUNZx+/IMk/TvIHsXfOXZsHYe9c+b/CS5Kzf97475LskrxzjPFvrnhK10JV/c3cfuqUJEdJfvVBX5uq+rUkr8vtTldqmWIAAACdSURBVN+fTfLTSf5Tkt9I8teT/EmSHxpjPHDfTH3O2rwutx+hjySfSvJjX/2ejQdNVf2DJP89yUeTrGcv/1Ruf6/PA71/7rA2b4z9k6r6O7n9TeK73H7w8BtjjH9dVd8Ye+e8tfmPeZ7vnWsRoAAA7ifX4Ut4AAD3FQEKAGCSAAUAMEmAAgCYJEABAEwSoAAAJglQAACTBCgAgEn/H20TlZAFG7mgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig= plt.figure(figsize=(10,10))\n",
    "plt.imshow(np.rot90(c_reshape[:,:,13]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "93750815",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from skimage import measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "73d8a437",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.037797786087383475, 0.14403888503599852)"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_reshape.min(), c_reshape.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "184b510f",
   "metadata": {},
   "outputs": [],
   "source": [
    "verts, faces, normals, values = measure.marching_cubes(c_reshape,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "57aa4ca6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2033, 3), (3990, 3))"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "verts.shape, faces.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "82ed7fa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f6feb8eac110467cbebeabb891be5954",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(14.060592…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<meshplot.Viewer.Viewer at 0x7fb17a4941d0>"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(verts, faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dc53937",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}