M2_SETI/D3/TP/TP_SETI_Kmeans/TP1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP1 KMEANS\n",
    "\n",
    "On nous propose de coder l'algorithme des kmeans afin de faire du clustering sur 2 classes puis plus de 2 classes.\n",
    "Plus tard, on utilisera notre algorithme pour segmenter une image sur l'information de couleur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy.spatial\n",
    "from skimage import io"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = 2\n",
    "dim = 2\n",
    "nb = 50\n",
    "K = clusters\n",
    "mean = np.random.randint(5, size=clusters)*2\n",
    "mean = mean.T * np.random.random(size=clusters)\n",
    "sd = np.random.random(size=clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fonctions à utiliser pour le clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_points(mean=1,sd=0.5, nb=100, dim=2, clusters=2):\n",
    "    \"\"\" Generates data\n",
    "    dim: dimension\n",
    "    nb: number of points\n",
    "    clusters: number of clusters\n",
    "    mean: mean\n",
    "    sd: standard deviation    \n",
    "    \"\"\"\n",
    "    size = []\n",
    "    # for i in range(0,dim):\n",
    "    size.append(nb)\n",
    "    size.append(dim)\n",
    "    points = np.random.normal(mean[0],sd[0],size=size)\n",
    "    for i in range(1,clusters):\n",
    "        points = np.concatenate((points,np.random.normal(mean[i],sd[i],size=size)),axis=0)\n",
    "    \n",
    "    return points, mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance(points,Pc):   \n",
    "    \"\"\" Returns spatial distance between two matrix\n",
    "    \"\"\"\n",
    "    return scipy.spatial.distance.cdist(points[:,:], Pc[:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kmeans(points = [0,0], K = 1):\n",
    "    \"\"\" Create K clusters from points\n",
    "    \"\"\"\n",
    "    # Initialisation K prototypes\n",
    "    dim = points.shape[1]\n",
    "    N = points.shape[0]\n",
    "    iter = 0\n",
    "    eps = 0.1\n",
    "    Pc_index = []\n",
    "    Pc_save = np.zeros([K,dim])\n",
    "    clusters = []\n",
    "\n",
    "    for i in range(0,K):\n",
    "        Pc_index.append(np.random.randint(0,N))\n",
    "    Pc = points[Pc_index,:]\n",
    "\n",
    "    while (np.mean(distance(Pc,Pc_save)) > eps and iter < 3):\n",
    "        iter += 1\n",
    "        Pc_save = Pc\n",
    "        # print(Pc)\n",
    "        # print(points[:,:Pc.shape[0]])\n",
    "        dist = distance(points=points[:,:Pc.shape[1]],Pc=Pc)\n",
    "        clust = np.argmin(dist, axis=1)\n",
    "        clust = np.expand_dims(clust, axis=0)\n",
    "        points = np.append(points[:,:Pc.shape[1]], clust.T, axis=1)\n",
    "        # print(points)\n",
    "        Pc = np.zeros([K,dim])\n",
    "        index = np.array([])\n",
    "\n",
    "        for n in range(0,N):\n",
    "            for k in range(0,K):\n",
    "                index = np.append(index, (clust==k).sum())\n",
    "                if points[n,-1] == k:\n",
    "                    # print(points)\n",
    "                    # print(Pc)\n",
    "                    Pc[k,:] = np.add(Pc[k,:], points[n,:-1])\n",
    "\n",
    "        for k in range(0,K):\n",
    "            Pc[k,:] = np.divide(Pc[k,:],index[k])\n",
    "\n",
    "        # print(Pc)\n",
    "        indice = points[:,-1]\n",
    "        points = points[:,:-1]\n",
    "    return Pc, indice, points\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors=['red', 'green','yellow','blue','purple', 'orange']\n",
    "def visualisation(points, index, Pc=[0,0], K=1):\n",
    "    \"\"\"Visualisation function of a dataset and its K clusters\n",
    "    \"\"\"\n",
    "    if(points.shape[1]==2):\n",
    "        # for k in range(0,K):\n",
    "        for n in range(0,len(points)):\n",
    "            plt.plot(points[n,0], points[n,1], 'o', color=colors[int(index[n])])\n",
    "        plt.plot(Pc[:,0],Pc[:,1],'ko')\n",
    "        plt.grid(True)\n",
    "        plt.axis([min(mean)-1,max(mean)+1,min(mean)-1,max(mean)+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def img_2_mat(my_img):\n",
    "    \"\"\" Reshaping 3D img NxMx3 to 2D matrix N*Mx3\n",
    "    \"\"\"\n",
    "    mat = my_img.reshape(my_img.shape[0]*my_img.shape[1],my_img.shape[2])\n",
    "    return mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mat_2_img(mat,my_img):\n",
    "    \"\"\" Reshaping 2D matrix N*Mx3 to 3D img NxMx3 \n",
    "    \"\"\"\n",
    "    img_seg = mat.reshape(my_img.shape[0], my_img.shape[1], my_img.shape[2])\n",
    "    return img_seg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kmeans_image(path_image, K):\n",
    "    \"\"\" Clustering an image and changing pixels to its closest cluster\n",
    "    \"\"\"\n",
    "    my_img = io.imread(path_image)\n",
    "    imgplot = plt.imshow(my_img)\n",
    "    Mat = img_2_mat(my_img)\n",
    "    \n",
    "    Pc, index, clusters = kmeans(Mat, K)\n",
    "\n",
    "    for k in range(Mat.shape[0]):\n",
    "        Mat[k,:] = np.floor(Pc[index[k],:])\n",
    "\n",
    "    img_seg = mat_2_img(Mat, my_img)\n",
    "\n",
    "    io.imsave(path_image.split('.')[0] + \"_%d.jpg\" % K, img_seg)\n",
    "    imgplot = plt.imshow(img_seg)\n",
    "    return Pc, index, img_seg\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exemple\n",
    "### Clusterisation 2D\n",
    "On fait la clusterisation d'un exemple simple avec 2 nuages de points éloignés"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnhElEQVR4nO3df3Ac9X3/8efbBitR3BhjS6pj0Bk3OIWWb0nxgNvMNBg7iclMMf1+XSq+gpo0HpVQ2jQJnZBRm1AaTUm+TJ30m6TEoxJopS+Oy3eSOA0pX2zkdObbOMX0CzaQMRiDHCiRUgPuKAaB7ff3j9tzVqfd+7l3u6d7PWZu7vaze3dvraR9735+rbk7IiLSvualHYCIiKRLiUBEpM0pEYiItDklAhGRNqdEICLS5s5IO4BanHXWWf7Od74z7TAi/fSnP+Vtb3tb2mFESiK2l197mfFXxznlp06XzbN55M7KcfZbz041tkZRbLXJcmyQ7fgaFdujjz76H+7eNWuFu7fcY9WqVZ5VY2NjaYcQK6nYRvaPeG5rzu0289zWnI/sH6n7M9thvzWCYqtdluNrVGzAPo84prbkFYGkq/+ifvov6k87DBFJiNoIRETanBKBiEibUyIQEWlzSgQiIm1OiUBEpM0pEYiItDklAhGRNqdEICLS5hJJBGZ2t5lNmtkTMevNzP7azA6Z2X4z+9XQus1m9kzw2JxEPCIiUrmkrgjuATaUWH8lcH7wGAD+BsDMzgY+A1wGXAp8xswWJxSTiIhUIJFE4O7/DLxcYpONwN8F013sBc4ys2XAB4CH3P1ld38FeIjSCUVERBLWrDaC5cCPQssvBGVx5SIi0iQtM+mcmQ2Qr1aiq6uLPXv2pBtQjKmpKcVWg6zEtmtiF8PPDTM5PUl3RzdbztvCmretyURsUbKy36JkOTbIdnxNjy1qStJaHsAK4ImYdV8Frg0tHwSWAdcCX43bLu6haahro9hKG9k/4p1Dnc5tnH50DnX64PbBtEOLlYX9FifLsblnO75mT0PdrKqhncDvBr2H1gDH3P0l4EHg/Wa2OGgkfn9QJtJ0g7sHOf7m8Rllx988zvBzwylFJNIciVQNmdl9wOXAUjN7gXxPoDMB3P0u4AHgg8Ah4DjwoWDdy2b2F8AjwUfd7u6lGp1FGubIsSOR5ZPTk02ORKS5EkkE7n5tmfUO/EHMuruBu5OIQ6QevYt6GT82Pqu8u6M7hWhEmkcji0UCQ+uG6Dyzc0ZZ55mdbDlvS0oRtafRA6Os+MIK5v35PFZ8YQWjB0bTDmnOa5leQyKNVrj95uDuQY4cO0Lvol6G1g2x/Kh6NDfL6IFRBr49cLqtZvzYOAPfHgDQ7VEbSIlAJCTqfsxZ7WI4F8U12A/uHlQiaCBVDYlIZsQ12MeVSzKUCEQkM3oX9VZVLslQIhCRzIhrsB9aN5RSRO1BiUBEMqP/on62/eY2cotyGEZuUY5tv7lN7QMNpsZiEcmUqAZ7aSxdEYiItDklAhGRNqdEICLS5pQIRETanBKBiEibUyIQEWlzSgQiGaFZNyUtGkcgkgGadVPSpCsCkQwoNeumSKMlkgjMbIOZHTSzQ2Z2a8T6rWb2WPB42sxeDa07GVq3M4l4RFqNZt2Upujk7KjiuquGzGw+8GXgfcALwCNmttPdnyps4+4fC23/h8C7Qx/xmrtfXG8cIq0s7jaZmnVTEvVzRN5lKYkrgkuBQ+5+2N3fALYDG0tsfy1wXwLfKzJnaNZNaYp5LIgqtvx95WtnZpuADe6+JVi+HrjM3W+O2DYH7AXOcfeTQdkJ4DHgBHCHu38z5nsGgAGArq6uS3bs2FFX3I0yNTXFwoUL0w4jUtZi2zWxi+HnhpmcnmTpgqUMrBxgfc/6tMOapVn7Lbw/uju62XLelrL7I2u/07AsxwbZjq9Rsa3tW/uG/9g7Zq1w97oewCZgOLR8PfClmG0/CfzPorLlwfNK4HngF8p956pVqzyrxsbG0g4hVpZiG9k/4p1Dnc5tnH50DnX6yP6RtEObJUv7rVi7xDayf8RzW3Nut5nntuYS+Ttpl30XRieHPeKYmkTV0IvAuaHlc4KyKH0UVQu5+4vB82FgDzPbD2SOUi8ZqVSha+34sXEcP921VuMsanCcl6OKk0gEjwDnm9l5ZraA/MF+Vu8fM/tFYDHw/VDZYjPrCF4vBd4DPFX8Xpl71EtGKqWThsarOxG4+wngZuBB4IfADnd/0sxuN7OrQpv2AdvdZzRKXADsM7PHgTHybQRKBG1A96aVSumkofESGVns7g8ADxSVfbpo+baI9/0LcFESMUhrGVo3NGMkLaiXjERT19rG08hiSUXxvWl7Onp0b1qJVGvXWs3dVDnNNSSpCd+bds+ePVx+0eXpBiSZVPgbGdw9yJFjR+hd1MvQuqGSJw2au6k6SgQiknnV3tC+VAOzEsFsqhoSkTlHDczVUSIQkTlHvdKqo0QgInOO5m6qTlskgtEDoyz9/FLszw37c2Pp55eqB4HIHFbcKy23KKdeaSXM+cbi0QOjfOibH+LNU2+eLjv62lF+71u/B6gHgbSG0QOjVfWakeobmNvZnL8iGNw9OCMJFLxx8g0NUZeWoLl2pNHmfCIo1UtAPQikFWiuHWm0OZ8ISvUSUA8CaQXqCimNNucTwdC6Ic6cd+as8gXzF6gHgZQ1emCUvr19qU5ToK6Q0mhzPhH0X9TP167+GkveuuR02ZK3LuHujXerIUlmKJ6b5qbv3MTAtweYmJ4oWTff6Dlt1BVSGm3O9xoC9R6Q8qLmprlr3104M2/lWjxNQTPmtCk1186ePXsS+Q5pb22RCETKiWqQLU4CBeG6+WbNaaOTGWmkOV81JFKJahpew3XzasiVuSCRRGBmG8zsoJkdMrNbI9bfYGY/MbPHgseW0LrNZvZM8NicRDwi1YpreDVsxnJx3bwacmUuqDsRmNl84MvAlcCFwLVmdmHEpl9394uDx3Dw3rOBzwCXAZcCnzGzxfXGJK0tjRuKxDXI3rj6Rno6emKnKVBDrswFSbQRXAoccvfDAGa2HdhIZTeh/wDwkLu/HLz3IWADcF8CcUkLSuuGIqUaZK952zVcfvnlVb9PpFXYzHvJ1/ABZpuADe6+JVi+HrjM3W8ObXMD8JfAT4CngY+5+4/M7BbgLe7+2WC7PwNec/c7I75nABgA6OrqumTHjh11xd0oU1NTLFy4MO0wIrVCbH17+5iYnpi1vqejh+1rtqcQWWvst1J2Texi+LlhJqcn6e7oZst5W1jfsz4TsaUpy/E1Kra1a9c+6u6rZ61w97oewCZgOLR8PfClom2WAB3B698HHg5e3wL8aWi7PwNuKfedq1at8qwaGxtLO4RYrRCb3WbObcx62G2WSlwjIyPe09PjZua5XM5HRkZSiSNOYb+N7B/x3Nac223mua05H9k/crq8c6hzxr7sHOo8vb4ZsWVVluNrVGzAPo84pibRWPwicG5o+ZygLJxsjrr7dLA4DFxS6XulvWSp8XV0dJSBgQEmJiZwd8bHxxkYGGB0NFuTvZWalE7zFEklkkgEjwDnm9l5ZrYA6AN2hjcws2WhxauAHwavHwTeb2aLg0bi9wdl0qay1Pg6ODjI8eNFB9HjxxkczNZBtNTBXt1bpRJ1JwJ3PwHcTP4A/kNgh7s/aWa3m9lVwWZ/ZGZPmtnjwB8BNwTvfRn4C/LJ5BHg9qBM2lSWbihy5EjMQTSmPC2lDvZZusKS7EpkZLG7PwA8UFT26dDrTwGfinnv3cDdScQhc0NWRtH29vYyPj4eWZ4lvYt6GT8WEWfQgyncCwvUvVVm08hikRhDQ0N0dhZVU3V2MjQUfxDN0hiIQjfWcldYacQs2aK5hkRi9PfnD5af+MQnmJycpLe3l6GhodPlxbI4BqKwPu7704pZskWJQKSE/v5+li9fHjugLKxZE9BFqbU6LS7mj373o0oEbURVQ9LWkqwWacUeOnGxHX3tqKqI2ogSgbStpG8K34o9dErFprEG7UOJQNpW0oOtsjQGolKlYsvylYwkS4lA2lbSVTmt2EOn/6L+GbdxDcvylYw
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "points1, mean1 = gen_points(mean,sd,nb,dim,clusters)\n",
    "Pc1, index1, clusters1 = kmeans(points1,K=K)\n",
    "visualisation(clusters1, index1, Pc1, K=K)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avec un K ne correspondant pas au nombre réel de groupes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApuUlEQVR4nO3df3Dcd33n8efbxlYQPhzbsYxxYskhcZq0hlAyIT1mih27ENohTmkKTpXUUDQqbcNBqDuE0RzQNJoLPaemLfSoRg1JsQ7jpgWcNkwucmw6UzDFuQtxEs6JY1shbpDATswJJfKv9/2x31W+Wn2/q939fnf3u7uvx8yOdj/fH/v2V/L38/38NndHRERa15x6ByAiIvWljEBEpMUpIxARaXHKCEREWpwyAhGRFveaegdQifPPP98vueSSeocR6ec//zmve93r6h1GJMVWGcVWmXrEduLlExz72TFOnT3F/LnzWfH6FSx+7eLMxFeqasX26KOP/tTdl87Y4O4N91q9erVn1Z49e+odQizFVhnFVplax7b98e3e3t/ufJapV3t/u29/fHsm4itHtWID9nvEPVVVQyLSFPp29zFxemJa2sTpCT72rY/R9fku5vzpHLo+38XQgaE6RZhdDVk1JCJS6LmTz0WmH3/5OMdfPg7AyMkReh/oBWAFK2oWW9apRCAiTWHlwpUl7TdxeoK+3X1VjqaxKCMQkabQv76f9nntJe0bV3poVcoIRKQpdK/pZuC9A3Qu7MQwOhd2suS1SyL3LbX00CrURiAiTaN7TTfda7qnPg8dGKL3gd5pjcjt89rpX98Px+sRYTapRCAiTSuqlDDw3oFpmYWoRCAiTa6wlCAzqUQgItLiUskIzOweMxszsyditpuZ/ZWZHTKzx83sl0PbNpvZM8FrcxrxiEjtDR0YmjZwa3h0uN4hSYnSKhHcC1xXZPt7gEuDVy/wPwDMbDHwGeDtwNXAZ8xsUUoxiUiN5BtlR06O4DgjJ0fY+vRWjeJtEKlkBO7+r8CJIrtsBP4+mO5iH3C+mS0H3g087O4n3P1F4GGKZygikkFR0ztMnpvUwK0GUas2ghXAj0Kfnw/S4tJFpIHEDdDSwK3G0DC9hsysl1y1EkuXLmXv3r31DSjG+Pi4YqtAubENjw4zeGSQsckxOto66FnVw4ZlGzIRWy1lJbaOtg5GJ0cj07MQX5SsXLsotY6tVhnBMeCi0OcLg7RjwNqC9L1RJ3D3AWAA4LLLLvO1a9dG7VZ3e/fuRbGVr5zYhg4Mse0726aqIkYnR9n27DYuv+LyqnQTbJbrVk13L7l7xsCttjlt3P0bd7N2zdr6BVZEVq5dlFrHVquqoV3A7wa9h64BTrr7C8BDwLvMbFHQSPyuIE0kVtx0w6qPrp+ogVtbVm9R//0GkUqJwMy+Su7J/gIze55cT6B5AO7+JeBB4NeBQ8AE8KFg2wkz+zPg+8Gp7nD3Yo3OIqqPzqjCgVtZrXaRmVLJCNz9plm2O/BHMdvuAe5JIw5pDSsXrmTk5EhkuoiUTyOLpeFETTc8NZGYiJRNGUGJCkdNaqBM/WgiMZF0NUz30XoqnMo2vNydbj71oYnERNKjEkEJ1EtFpLGpRF+cMoISqJeKSOOKmgep94FeTYoXooygBHG9UdRLRST74kr0g0cG6xRR9igjKIF6qYg0rriS+9jkWI0jyS5lBCVQLxWRxhVXcu9o66hxJNmlXkMlUi8VkcbUv74/cgH7nlU9dYwqW5QRiEhTyz/A9e3u47mTz7Fy4Ur61/ez4rhmvM9TRiAiTS+qRK+5kF6lNgJpSuo3LlI6lQik6WgkuEh5VCKQpqOR4CLlUUYgTUcjwUXKo4xAmo5GgouURxmBNB2NBBcpTyoZgZldZ2YHzeyQmd0esX2bmT0WvJ42s5dC286Gtu1KIx5pbRoJLlKexL2GzGwu8EXg14Dnge+b2S53fyq/j7vfFtr/o8BbQ6d42d2vTBqHSJhGgouULo0SwdXAIXc/7O6ngB3AxiL73wR8NYXvFRGRFFhuXfkEJzC7EbjO3XuCz7cAb3f3WyP27QT2ARe6+9kg7QzwGHAGuMvdvxHzPb1AL8DSpUvftnPnzkRxV8v4+DgLFiyodxiRFFtlFFtlshwbZDu+asW2bt26R939qhkb3D3RC7gRGAx9vgX4Qsy+nwT+uiBtRfDzYuAo8KbZvnP16tWeVXv27Kl3CLEUW2WyFNv2x7d757ZOt8+ad27r9L4dffUOKVaWrluULMdXrdiA/R5xT02jaugYcFHo84VBWpRNFFQLufux4OdhYC/T2w9EJBC10tbWp7dq+gxJLI2M4PvApWa2yszmk7vZz+j9Y2a/ACwCvhtKW2RmbcH7C4B3AE8VHisi0SOmJ89NasS0JJY4I3D3M8CtwEPAD4Gd7v6kmd1hZteHdt0E7AiKJ3mXA/vN7AfAHnJtBHXKCIaALnKXpCv4LJIdGjEtkYaGoKsL5szJ/Rwq/96VyqRz7v4g8GBB2qcLPn824rjvAGvSiCGZIXLt0PmnrZHgM4C6IEo2rFy4kpGTI5HpWTB0YGjanP83L7+Ztaytd1jNbWgIenthIrh3jYzkPgN0l37v0shiAPp4NRPImwjSRbIhasR025y2TIyYVvtFnfT1vZoJ5E1M5NLLoIwAgLiitYrckh1RI6a3rN6SiYFzar+ok+di7lFx6TG0HgEAK8lVB0Wli2RH4YjprKyypfaLOlm5MlcdFJVeBpUIAOgH2gvS2oN0EZmNZnytk/5+aC+4d7W359LLoIwAyDUIDwCdgAU/B1BDsUhpstx+0dS6u2FgADo7wSz3c2CgrIZiUNVQSDe68YtUJl9dVdhrKAvtF02vu7vsG38hZQQikoqstl/I7FQ1JFIFQweG6Pp8F3P+dA5dn+9SN0rJNJUIRFKW71Of7045cnKE3gdyg3xUVSJZpBKBSMqi+tRPnJ5Qn3rJLGUEIilTn3ppNMoIRFKmPvXSaJQRiKQsqk99+7x29amXzFJGIJKyqDmBBt47UPeGYvVkkjgtnBFo/QFJX/5me8s/3QLAV973FY5+/GgmMoHC2UF7H+hVZiBAS2UE4Rv/BcCHyE0057y6/oD+U0jlsnyzVU8mKSaVjMDMrjOzg2Z2yMxuj9j+QTP7iZk9Frx6Qts2m9kzwWtzGvHMlF94Jn/jPw6cLthH6w9IMlm+2darJ5OqoxpD4ozAzOYCXwTeA1wB3GRmV0Ts+jV3vzJ4DQbHLgY+A7wduBr4jJktShrTTFELz0RR9z6pXLVvtkluqvXoyTQ8OpzZEpJMl0aJ4GrgkLsfdvdTwA5gY4nHvht42N1PuPuLwMPAdSnEVKDU/4jq3ieVq+bNNq7aaXh0uKTj69GTafDIYGZLSDJdGlNMrAB+FPr8PLkn/EK/ZWa/CjwN3ObuP4o5dkXUl5hZL8FCwkuXLi1rQqtrrungvPNGi+5z9mwbBw/ezNhY6eeNMj4+ntnJthTb7IZHhxk8MsjY5BgdbR30rOrhmtddU1JsNy+/ma3/byuT5yan0trmtHHz8psT/9v+eN8fR95UBw4PsGHvhlmPX8EKbnvTbTP+bSuOr6jadR+bHItMf+7kc5n4XWflby5KrWMzd092ArMbgevcvSf4fAvwdne/NbTPEmDc3SfN7PeBD7j7tWa2BTjP3e8M9vuvwMvuvrXYd1522WV+8ODBMqIsXJweYD7wn4AT5EoC/aQxDfXevXtZu3Zt4vNUg2IrrnCOIMg9Nd/2ptu48wN3lnyO8FTM/ev7U+kxNOdP5+DM/L9qGOc+cy7x+avhDXe9gdHJmQ9gnQs7Ofrxo7UPqEAW/ubiVCs2M3vU3a8qTE+jaugYcFHo84VB2hR3P+7u+cekQeBtpR6bjqiFZ+4BfgqcA46itQgkrrF38MggUFodffeabo5+/CjnPnMu1W6jcdVLHW0dqZy/GnpW9WhgXYNIIyP4PnCpma0ys/nAJmBXeAczWx76eD3ww+D9Q8C7zGxR0Ej8riCtCrrJ3fB145docY26Y5NjVesaWmoDcFwdf8+qnsj9s2DDsg2ZHFgnMyVuI3D
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "K = clusters+1\n",
    "points, mean = gen_points(mean,sd,nb,dim,clusters)\n",
    "Pc, index, clusters_ex = kmeans(points,K=K)\n",
    "visualisation(clusters_ex, index, Pc, K=K)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et un exemple un peu plus complexe avec des intensités différentes et relativement proche"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = 5\n",
    "dim = 2\n",
    "nb = 50\n",
    "K = clusters\n",
    "mean = np.random.randint(5, size=clusters)\n",
    "mean = mean.T * np.random.random(size=clusters)\n",
    "sd = np.random.random(size=clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjaklEQVR4nO3dbWxc13kn8P/DEUVprJVcUyKbtc1hDMcEjGibNEZrw7tYyVIBR4AadHebOhgL3t0IA9lNIQMVugUG2NTAzi5QuK30IXZCKE4FcVBDQNI2yioIIkVMNkW9jV04pl+q1FA5XHcDEaJb2jQtmi+nH+5ccl7uue8z9557/z+AoObtzuHV8LmHz3nOOaKUAhERmWsg6QYQEVE0DORERIZjICciMhwDORGR4RjIiYgMty2JN7399tvVvffem8RbG+ODDz7AbbfdlnQzUo3nyBvPkTeTztErr7xyUym1r/P+RAL56OgoXn755STe2hjT09M4cOBA0s1INZ4jbzxH3kw6RyLScLqfqRUiIsMxkBMRGY6BnIjIcAzkRESGYyAnIjIcAzkRkeEYyImIDMdATkRkOAZyIiLDMZATERmOgZyIyHAM5EREhmMgJyIyHAM5EZHhGMiJiAzHQE5EZDgGciIiwzGQExEZjoGciMhwDORERIZjICciMhwDORGR4RjIiYgMFzmQi8gOEfkbEfmpiLwhIs/E0TAiIvJnWwzHWAHwiFJqSUQGAfxYRL6rlHophmMTEZGHyIFcKaUALDVvDja/VNTjEhGRP2LF4YgHESkAeAXAvQC+opT6bw7PqQCoAMC+ffs+c+HChcjvm2VLS0vYtWtX0s1INZ4jbzxH3kw6RwcPHnxFKfVA5/2xBPLNg4ncDuDPAfyOUup13fMmJibUtWvXYnvfLJqensaBAweSbkaq8Rx54znyZtI5EhHHQB5r1YpS6p8BXAXwaJzHJSIivTiqVvY1e+IQkZ0Afg3A30U9LhER+RNHj/xjAK6KyGsAfgLg+0qp78RwXCLKuHodGB8HBgas7/V60i0yUxxVK68B+HQMbSGiHKnXgUoFWF62bjca1m0AKJeTa5eJOLOTiBJRrW4FcdvysnU/BcNATkSJmJtzvr/RYLolKAZyIkrE2Jjz/SJWMFdqK93CYO6OgZyIeu6pp4Bt26wgvW2bdbtWA4rF9ueJWAG8VWu6hYOjzhjIiSgytwD71FPA888D6+vW7fV16/Zf/RUwOQmUSlYAL5W6g7it0QD27gUef5y9dScM5EQUiV19oguwk5POr3v+eaunfeSIlWaZmwMKBf37LCx032f31vPeU2cgJ6JIdNUnJ09aQdXuiTtpNKyAbl8E3J7rdgy3C0keMJATkW9OPV9d9cnCghVUe61QYBljHOuRE1EOOE3gefxxK6jHuPZeYLpevO4Ck0XskRORL04pFADY2Oh/W2wDLhFMV96YRQzkRORLGnu4uotIsWiVN+YFAzkR+WJKD7dUsipl8rReCwM5UY4FKdtzmsCTNiLA7Gy+gjjAQE6UW171353KZaunK9LfdgYxMJCvskMbAznlQ95njDgIs/pgWnq627cDg4Pd96+v56+GHGAgpzwI2vXMCd3gpdegpi5XPjzc29RL61T+F14AvvEN55mgeashBxjIKQ+48LUjXUD2GtR0ypUXi8CZM8ATT7hPsw+rULBy3+fPW7ePHbP++1hDbmEgp+wL2/U0TNDskS4ge5Xt2bny1h6yvZ7KuXPhptl7mZhwXjRLl683pcImLgzklH1hu54G0WWPLl8e0b5GF5DLZe+LQrls9ZA3NraqRHQThuLw5pvOi2Yp1R3M81ZDDjCQUx6E7XoaRJc9Onv2HtfXOQXksEMK1h84dQDjsELLePN2bynlfDHKEwZyyj63rmdG6LJE8/NDgY+luyg88YR72qZYrAOoAGgAUM3vFfQ6mJdK3RejvGEgp3xw6npmiC5LNDKyEvhYuovC+rp7D/3DD6sAOnMrywCqbdfQ4eHATdLK2B9WoTGQE2WALnt0/Pj1wMfyM3TgVPSzsaEbPJ7DkSNb19B33w3cJEcimfvDKjQGcqIMKJeBhx5qv++hh4DDh+cDH+vee/09r3O3+4EB3RVgrG2XIL9jzIWC9yxSBnELAzlRBjz1FHDlSvt9V64Ap0/7jMpN9Trwgx/4f37rgKhSNQCdM4KKAGptJYl+1mwpFq1Sxo0NKx3jJENFR5ExkBNlgG5fzIsX/3Wg41Sr4TeJUKoMYBJACYA0v08CKG9OEqrXtwZTdWuJDw+3p0xyUHQUGXcIIsoA3SScjY1gK1xFnyNVbn61W1+3JvS8/z7w0Ud226z1UnbvtvLmY2NWcO5Ml9i3q1Wrfbrn5RkDOVEGFArOwXxgQMHqHfszNta7fTadJvSsrgK7dgE3b7q/tlxm4HbD1ApRBlQqzvcfPfr/N//tZwp/EmuOZ2ylhEQwkBNlwHPPAU8+ubVgVaFg3X766bcB+J+tWS63L3zVj7XHOWgZHQM5UUY89xywtmYF6rU167bNawFIu7cuAjz//FaaJuzAp18ctIwHAzlRDrgtANnaW4+iUAg2azODKyUkhoGcO8dQDujSFwMD1tKwUVcttOu+P/95f8+310dhEI9HvgM5d46hjOnsl9jL2OoGMaOsHX7bbd1rkF265O+1HOCMV74DOXeOoQxozW8fO9beL3n22QnU690LQEbZxWdgwBpI/drXrJ7+3Jz1K1Ov+w/QHOCMV74DeU52jqHs6sxvdw5OrqwUNvslrQtAbmyEf8+77wYeftj5j9k77vB+vdMAJzOc0eQ7kOdg5xjKNj+78jj1S6J8xO0euNMfs0B3Cmdw0BoE1S0FzwxndPkO5FzEgQzn549Hp6AdZeLPHXfo3/fdd7v38PjGN6yZm7ql4JnhjC5yIBeRu0Xkqoi8KSJviMjJOBrWFznYOYayzatnPTS07tgvsT/6cb/v2FjwPTyY4Ywujh75GoDfVUrdD+BBAL8tIvfHcNz+yPjOMZRtTj1rezZmqQScOnWtK41h56Kr1XC79SwsuP8xGzTfzQxndJEDuVLq50qpv23++30AbwG4M+pxicib0x+V589buebZ2faNJZxy0e+/b+WwgygU9H/MAsHz3cxwRicqxjm4IjIO4EcAPqmUeq/jsQqsnVixb9++z1y4cCG2982ipaUl7Nq1K+lmpBrPkbfWc/TYYw/ixo0dXc/Zvfsj7Ny5gRs37I2avRZYUbh69YeOj+jeY3T0Fl588SXtES9fHsHZs/dgfn4IIyMrOH78eqjdjcIw6XN08ODBV5RSD3TeH1sgF5FdAH4IoKaU+pbbcycmJtS1a9died+smp6exoEDB5JuRqrxHHlrPUcDA85rp4hslSPqntPKnpXpxM97pI1JnyMRcQzksVStiMgggG8CqHsFcSJKhp9ctFde2ivlwXx3MuKoWhEAXwfwllLqj6M3iYh6wU8u2mvw1Kuoi/nuZMTRI38YwDEAj4jIq82vIzEcl4hi5KfatnM98kIBOHFia/DUq6iLFb3JiKNq5cdKKVFK/Rul1KeaXz6XziGjcV61cbyqbet1axVDezGt9XXrdpD/Wlb09l++Z3ZSeJxXnUmcZWkmBnIKh7/xmcRZlmZiIKdw+BufSaw6MRMDOYXD33gjeQ1rsOrETOYHcg64JYO/8cbxM6zBqhMzmR3IOeCWHP7GG8fvsAarTsxjdiDngFuy+BtvFA5rZJfZgZyfTCLfOKyRXWYHcn4yiXzjsEZ2mR3ITfxkZm1w1oSfx4Q29gGHNbJrW9INiMT+BFarVjplbMwK4mn9ZNqDs3Ze3x6cBdLbZjcm/DwmtLGPyuVc/tiZZ3aPHDBrwM3EwVm33mw/f56wvWoTzzlRQOYHcpOYNjjrVd7Zr59H046Ry5e9X2vaOScKgYG8n0wbnPXqzfbr59G0456zZ71fa9o5JwqBgbyfTBuc9erNBv15wqZHNO0Ymvexp6Np55woBAbyfjKtbMBPb3bnzq1/Dw/rf54os3A17VgZGfF+rWnnnCgEBvJ+M2lwtlYDtm9vv2/7dut+OzAvLGw99uG
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "points, mean = gen_points(mean,sd,nb,dim,clusters)\n",
    "Pc, index, clusters_ex = kmeans(points,K=K)\n",
    "visualisation(clusters_ex, index, Pc, K=K)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exemple de clusterisation sur une image\n",
    "On souhaite pouvoir changer les pixels vers les le centre du cluster le plus proche.\n",
    "\n",
    "On observe ainsi pour un nombre différent de clusters :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x1a9ff9b4dc0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAD8CAYAAAAVOD3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOy9d7wkd3nm+/2FSt19wuQ80mg0GuWAIsoIBQQCgYgCY4yx2cXYXnv37q733g1317vrzfde413bBBEkhCRAQoByQDmO8ihN1uR4UodKv3D/qO5zzoxGoBESkpbzzKenu6v7VFVXV7/1hud9XuG9ZwpTmMIUpvBqyLd7B6YwhSlM4Z2KKQM5hSlMYQqvgSkDOYUpTGEKr4EpAzmFKUxhCq+BKQM5hSlMYQqvgSkDOYUpTGEKr4G3zEAKIT4ghHhZCLFGCPEXb9V2pjCFKUzhrYJ4K3iQQggFrAIuADYDjwOXe+9feNM3NoUpTGEKbxHeKg/yFGCN936d974ArgEufYu2NYUpTGEKbwn0W7TeBcCmSc83A6e+1ptnzpjmFy+a/xbtyhR+O+ERr/WKByEEHvDeIRAIUf2FEALrHFK81l/vvQ08CAHOOaTS4BweqoWvgqj+hjcStb2e/ZnCG8VTz7yw23s/a9/lb5WB/JUQQnwZ+DLAooXzePCOq9+uXZnC/5aoDOR+zZSUOGsxxqCTBGEteZ6jA43SmqIsCMIQayzeezyg5KvXJADpPViLD0Pa7ZRGFIJUCATWT+yF675f4Kpnr9fejdvSKQP5VqI2+/hX9rf8rTKQW4BFk54v7C4bh/f+68DXAd5z/FFTDeFT+I1BKIUtCqJaDVeWKKUIA413Fi0VpclBODSVp+m8R9j9r0sKidQKhyfUEinAeQtIhAe6nmmVyxKA+w19yim8GXirDOTjwDIhxBIqw/gZ4LNv0bamMIUDg7UIKfHGYKwlyzLysmRgWj+7t27h1ptuYvfu3QRBgJSSoijQ+tU/FRVoLv/s5wCBlJI4jkFoTFGghEIohZjkRSJ8N/z+DX/eKbxhvCUG0ntvhBB/DNwGKOAK7/3zb8W2pjCFA0We58SNBmmrRdhoUI6O8rMbf8LRhx/GyaeczOWXfxYdBHjnqqyh90j56nqmMSXf/8E1eCmwxnDCCScQRTFzFy5isL8P3/UkwXfznbLyKL2YMpLvErwlNJ8DxXuOP8pP5SCn8ObitXOQ3ntuuukmlh9+OEuXLq0CX+dQ3qNDjS0KlFLYskRJBUoh/H5CYw+lt+AdUmseeuB+Nm3axMEHH8xpp502XgjqZR894JA4IV6jiLPfj9HFlEV9K1GbffwT3vuT9l3+thVppjCFNwOia2iccwRBQJZlVU4xDCnzrAqlux6gtRalFEIIjjrqKAYHB8fNjlIKPFgncQR4LxFaYxF47xGo/Wzco6XAWoNAcOaZZ9Jut3HOcfXVV9NoNLjkwx8GD0pKhJRIpegUJUEQAFAUBVJKtNa8E5yVKeyNKQM5hXc1nHMIUdFzsjxHhyFFWZK3WjSSGOc9WZoSxzE6DDFliQCWLFmy13o84IXCI/Hj+cZulWX8nr0eC2+r50rj8EgktXoDgM987nOYosB7T7PVYrB/gCLLQVmE0qgkIW82UUoRNRq0hoeJouitPlxTOEBM9WJP4V2NngdpjEGHIWWeEwUBYa1GWZY470mSBK8UWadTGdP95BOhMn1WCNz4Te5zv/dj333skXhUFT53b1KFFNaj63Vu+PGPefjRRwiSCB2FaK3Jm83x7RbtNvVG4zdxuKZwgJgykFN4V8N7j3OOqF5nZM8eao0GRZ7jvEdHEQhR5RelJAiC1zSO3bV1b+513bxwFatRSJyQWCHxQuKpthnGCUVWcPa57yNMalz1ne/x4sursN6Nh9VKKbTWZGn61h6oKbwhTBnIKbyrIUSVIyzSlGnTplHkOSqOGd2zhx98//s0Wy2iRgNbFNUfOFe10uwXDol53bce/duP3ySOyrMUKCQCjGP+/IUce/QxHHXsMWxcvx4RhCAEeZ4ThiHNVou4VvtNHbIpHACmcpBT2Ad+P89+dQX1Ve/w+7ywr00S+1n2S1boX+Ml71xVYBECYwxRHLPhlVd44IEH+NQnP0lcr9MZG8M7RxzHVWdM10Duu3mJx++nWr03d9Hv8xd6vKlRjFerq/foMMaZEh2F+CBg+eFHYP1LPPXww8w96CDmTJvWzYl6TJFPfLqpWs07BlMe5BT2A4cQDkTVEucEOCFBBxipMELglMJJgcGDkDgvEUKDV0ihUF0PSnqBLw2BVGjnoSjRSjG+eifQMkCh0EKhvETpAIxFAVoIvLVIJSmsqXw2KRFIJIwTuLVSlHnG9m1beO7Zp7n0so+ilCRvd1BKEyUJRVninOvmDOl6fL2bwJYO5STSQiBDhJfgBVJqHALjPU4CgaLA4mTlN3rnEVJhHAgVYZyi8IrcgpWavHQUrTa1eoNT3nMiRbPJ1X/3dTauW4t3hkY9qTiToptT9Q7lBdJJFAr8FMXn7cKUgZzCqyCwVYXWezweJ6riRWYthbUYgEDhpQTVDTSFxDmBsRZbOrK8IAhivPWouI4zBuup+pSNxTlPECUoqSk6Bd4ATuCsB1O1/JmiRGmNklURJojibo+0q9r4qmZnyjRFeE+9ViOJI04++UQaSVyJUkgJQmCtQ0oNKLyvgumJPmmJQ6HiOjqsIZ0mTwvKwiJkQFaUVQgtIC8LOlmHMAqwziCEx+QZJi8BQWEdKooROsCiMF4gVIAKwvHPdvY55/L+s8/imquu5OUXVoISeG8xpqAocgIdooKoam+0ouuZTuHtwFSIPYX9QO1FcK5+nh4tqQyULZGmpMhSnDXESR3jLBJFUI/BGKI4ZKzTJtAa6RxKB9iyJGo06DSb6CgisyXOWcIkBMB5TxBq0jwl6avjWgXGGqwQlMZQiyOyNEULidaS9libWr2OVIpdO3ZQr9dpNProGxyk2W6TBPGr+Nh+nKYz6QVRfVqTp5TOEugAHSiCIKA0BUkUUNqCRqBodXL6ZkzHtJqEQuMDRRIGqHqd5miTOKmRZylSKuR+BC5QirKTcuo5ZxP11fj5T39Kf38fi5cswZS2qopbS14WBCpESoXbH0l9Cr8RTBnIKbwKDtm1Hz2+X5VZK/KMuFbDpB2kEtSlQwYK6wpkGJNnGeRVGKuDCBkIkGCcAaUhCOgUBaqWYLzFGoNUAiJN1moSKIVzEhUFpJ0OMgjITImQgloS48oC6RxBHGBKQ9LXYM+uXfzkJzdQGssXfv/3Kq/QQ61Wx5t9DMtkm8hEHlL67qOw8jbLskBLDcLQHt2D1pJVK58j8o4gDHDTp7Nl0yao1Vmy/BiiWj/tPTup9U/D2wKNQwrZlVLrqvd0t++FgChAIDjqmGN45tmn2fzKJhbOXwDGEUQJxriux5xgvZ1KSb6NmDKQU9gL4w163cqEECC8R1T9JWgcQnpkpPGtNlIGuNIipEPakjBMMM6gVIDJc4K4jyzNKIwhDGOMMVWI6w1BpMA60rRJUk8waQdnHfG0GZTDw3gBYRRTlAW2NCgBcbcV0DuPcZbHHnuUI448ksMOP5wgCCmNgS55vNJ03G+z4fij3js8HiEFrbQimAtfYjotEmH4yTU/5KWnnqA/jgiBhfMXsGndOqwK+Ed/9d/AG/r6aqRjQyR9/Ujv8d50U4d+0lYFpfc4KTF4pNacdvKp/PAHP8CVltPPOYcyTVFRTNJXJytylFbd4z+FtwNTOcgp7IMeXaUiPnsvurqHkIQhPu3gs5R853au+dbf843/+p958v670aIgH9mOy5uIfAyRjVGThtbOLdxzy08JsZC1iAOJzVpEEnyZU6YtYulpD+8mCjRJFLB93WqSOCKUikBIhIOyKFB48k4HiSMMFa3RYbbv2MYJp5zCjFmzqvcoRRDHVSHHy3H75Pe6Vf8LfGV8vEd6j8sz+pIIV2Y888gD/M1/+0/c9MMf8OKKh+nTHtcew4wNs33dKgJXIMuMm6+7BlF2uPHaqyh9SXt4FwqLxqK8RXbX3xM6M95jBMgwwBYly44
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAD8CAYAAAAVOD3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9Z4xlaZrfif1ec8y14SMy0pfJ8rbLdLUf1+RwDFdLkcNZUhIhCuIHgYIAQZAI6IMAAQIoCPogAdIuZ3c4u8ulaIYc0+NnuqfNlPeVmVXpXWRk+IgbEdcd8xp9eM+NjKyqHnZPV3XWbMc/cfPeuO6ce8xzHvN//o/w3nOAAxzgAAf4KOTdXoEDHOAAB/is4sBAHuAABzjA98GBgTzAAQ5wgO+DAwN5gAMc4ADfBwcG8gAHOMABvg8ODOQBDnCAA3wffGoGUgjx80KIC0KIy0KIf/ZpLecABzjAAT4tiE+DBymEUMBF4OvAIvAG8F947z/4xBd2gAMc4ACfEj4tD/J54LL3/qr3vgD+LfCffUrLOsABDnCATwX6U/reI8DNfX8vAp//fm+enp72J0+e+JRW5QD/o4EXP8Cb3OjNIBwgAIX3AgT46jWPQyLC316AFwhZfUwweiPeewTha7wHIT60DqP3CfDCVcv3gEAgEYjqxX2f+UF+Bv4HfN8BPgm89dbbG977mQ8//2kZyP8khBD/BPgnAMePH+f11169W6tygL8GCIZMBkOGQHr2GUwfbtICFoTFA0ZKPBLhgykUeHx1b5xBSoVzIGREaQxKCBQOpSVYgxQeZ0uE9HjncNajdIRHIaSkLC3KC6SSoCK8FxgPTjhQDrxDOJBOooRGCI3we79mH8QdD8NPC+saPuCDoT3Apwal4xsf9/ynZSBvAcf2/X20em4P3vtfA34N4NlnnzloCD/AX4rg65nKvQOQe3bFCw8+eHCies0R4azCS4EUBsUQSVl5cwnKKrwRWJ8jkwKt+wgKfNHD9IZ0t7dwtsDlGTpS9Lq79PsDpmfnEDrGe5iam0S1NbgU7BSYCZSto6RExh4vLU4VoDwOg3MGCUgpKkNJ+A2V0d/7XdUt2MbKExYHp8jdwKdlIN8ATgkh7iEYxl8F/sGntKwD/I8eHnAI4faeceLOV4PhFCGURiG9IraA83gFXihKGQyPEglCefA5stjGrt1ge+kMrlwj662SdXfpbG2gJURSUEtj8mxAUZTcuhqhlEZIxVqSMHPyAWQ0SWvifqKJh5FqDkEbsgikwusUF5U4UeJUgUYCCrFXHA2htIDgIaPBa0TlLQP4/TH/AX6s+FQMpPfeCCH+KfAngAL+pff+/U9jWQf4yYDwVF6WxwsbvEZ8ZUNkFX5rELoyNB5BCd6BFaAihAvepil3KbtrlL1lejfP0Ll5lmxnGaUMpSuYmphgrj2OwLHb2SIbZkggEYKJRhsEdLtdysJx7f1L5BZqtfNMTLzD2PgJmtOnSGcfRESTaDSlAascQgbjiNd4HAiHx1Wh/8ghdohR3pMDs3i38anlIL33fwj84af1/Qf4CYIX4BXBIPrKaNiq9uErD0sivNrzvDyAGgZj4xOEFQgMyD7dtfe5fu5Fss51GG4TmZxmmtKYuJdSNhkbn0DgkXhEuo33VT5RCNJGA6UUupWRlYaalJR2iPN9tjtXWF05T7N9mmP3PUdr7iHisWNEURPhFF5GCB/tGXAweGHw2NuZA08Iqb3jNsnkwEzeLdy1Is0BDvDDQDAKOf1e3g5G+UeJrCrV3oEK9pJCpkhKtDO43irb518n37nEcOs8KlsiMj0Kr5HNWdqz91JrH6JfQiZVKMo4i56ZQSCwzqKlYFAalJLotI4yBi0VdRySkomxHsb2ENKyfOO7XL78IlPzj3Lk4a+Qtu7D5U2wEhFJnPegYpSKcFicK5GCfblIA0Iiqt98YCLvDg4M5AH+eqCi4oAAF+Glx4+8Lecx1qGlRGnAeIq8pEw8FF1sscHyW39KtvA2dbOKGK7Qii2+1qAvE4bCMJQ53g9wKOIoDbURNA7BYDAkTWsYKSjxGDxWeQwCISTSS6RX6FiRiggpBnS2V/FlxuLlbfqdTR7+3C8QTTwEqgleoISgKDzGO+JEo5TCeVMVo/xtL1JAyFIdVLHvBg4M5AH+GkHcvq/oMsJJBA6JRfkSYQDvSGSJKrZZu/ke27dOM1g9jchvomTO1Mw03e6QbAAu9thyi96mJZGOWM/gczDWIXSMTlK01Djr8F4GbqMQe5xIh8SLQO3xXiCJUNSZnTuFNYbV5XW2V2+wdO7bzJxcJR5/hCg9DlqTxIrI62ASnQCiqmptGNGVwm+WHBjIu4MDA3mAvx4QIwJ2MBbCEyrXPnhcStmQYzRDKId432P73HfprZxne/0aaa1EtBrUx48T1SfwskNqLc2moj9YYbe3yI4DlwxpTsxQa7bITEnpwSQCV3lzodACwThD5CTg8MJgpQlhvq8hZQuPY3KySa3WYWPtA3q9S8w+NGD2RB1cArIOXuFcBFIiR4lIrxlV7g/yj3cXBwbyAH8N4Ake1cibGh22VV3bGZwdIBMLbof+5g0GaxfonP1DbLZFO46J4wla04dJmpM4kdDSE0QuQ2RbDLcN9cLQNx12sAgKlJgBmaBkExnFgQAOSClxzmGsQSHQXgK2oqBLLApLhPAaFUGzVqM9punu9llavgLXXwXrmTr8ACqZR5Q1cDWUquECX6n6aSM+pOMAdw8HBvIAd+DD5YCPdHD4Ee/wzifFvk98nM8jqs/e4REJv++Tt8u44kNvCy+P2vhUSEd6VVFjPFI4hDSQbbF98wNuXD2D2b7ChF8nSRwuaiBdjeF2yaC3hY4FuCE+65Ltdii6hlRN0fARCRa/vkZ/c5uJucM0xjTb/W1qSR0nFHGS4BFYY5FSIClwaKyIsSQ4IfHC4USGkwXODdAuR8c14mic5cWrrGxucm9vmeMPvoBqnkS7ZF+D4qh3UbFnHA+i67uGAwN5gI/ACxeC2YouI3wwky4k3rAenAehBM57lHeoqgLrnEcqTVGaQKgGhPB47wKX0dvwNxbnDEIlWDSFKUmiOFBr9ooVVVHGC7AKryIsfYwoQKVYG5F6kA48GYOFt1j64Ltk2zdJxADnCuK0hrMFeXedXrZMUqvhRImgRHqD9JDGEcJJrHNICQKLsUPynXXMsE9vWDI3fxxUhLICpMJW1XNfFdWdcCBLkAKkxVMgtcXlFiEU9XSCyTFPVi6Q56tsXP4OabbL3GO/iGy2sc5gtCcSDuUFlBKEAhXoTQe4OzgwkAf4EPxtuQfhA83EB+EFKYJ/J2XoES69RUiP9CCxeOswxqFljdJadBRXBQgP1iFxIW0IgEVIhXEerxVOSTJj0djgxbnwuncSZISwCpB4LfDSYsixDuoywmd9hjffYf3Cd2H7fcboI5xFEJNnoe4syWjFAuyQ0Rrc9nVDK4uViihpMDU1RV4UbG1t4vvbICM6awadNIjrbeKkgYySsP4lIB1R4hHKYnyJQ6CVRKCABGsMRkgatRb3HDvO2voFGmaH7OZ7rAxT5p6fQbYP42SJ9QXSKaSvIWRc9YpbDozk3cGBgTzAh1AJQYyoJSNqTXV+FsMhSoOKPM5lJJECbysitURHCuFL8AWeBITES5BC473Dy6o4a2OkVuDAG4hUaBGM9uLrUQVXYS2oGHxRIKMEiQdnqCuHL7fYWXqTG+/8MfQv00gc3ggQUSXxMEoS+O8T/IfXvABsiaIg0Q5flkiThS2iIspejhnuMuxu4oXGBl8XHUeoKKUxMU3UmkT7CE+EcJrSlERIsmKXjc4yU2NNGvWYWDXZWtsikh3I38adn+Doc3+bmp9GiAbOGbwKBtcbEzigB3H2XcGBgTzAnai8xT1qSSXVNcqPpWkEZGD6SAqKrQ55PkDoCGPBeklzfJpIxdhiBxWlWCtQOsZ6h7MgvUIhwYUleQc4j/AOgcELhxQSLyQ2gmEpaOigqIOPkFZQExLldtm68SrL5/8cm92grvKKFqlRUlVf/INBAJGWlNmAlVs3wHu0sGilKMo+EoVwOnAwhUI
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAD8CAYAAAAVOD3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9WZBk+XXeCf7+y918jz0j96wdhZ0oUFxFiiJFiW3dMvWMNKNZrB/aTE/zNC+jp3kYszHrl34bs56WbHq6JTNRItUiKYgkSAIkCAprFYAq1F6VVZVb7BEevt31v83D9cjMKqA43UIVijTEl+YRHjfcw6973vvd8z/nO98RIQTOcY5znOMcPwj5Ue/AOc5xjnP8VcU5QZ7jHOc4x/vgnCDPcY5znON9cE6Q5zjHOc7xPjgnyHOc4xzneB+cE+Q5znGOc7wPPjSCFEL8XSHE60KIm0KIf/phvc45znGOc3xYEB+GDlIIoYA3gF8D7gHPAv84hPDKB/5i5zjHOc7xIeHDiiB/GrgZQng7hNAA/xr4+x/Sa53jHOc4x4cC/SH93UvA3Yd+vgf8jfd78Pr6erh+/dqHtCvn+GuPcP8LIN7z/aGHifDQ48J7HvNgpSQQhADi4fjgvX/+bINo74ez5xHa5wqBQEBg+XP7/ezxCA+E5Z+TEORyl0J7a3f4h76P9vk/dKfO8SHhO9/57nEIYeO92z8sgvz/CyHEPwH+CcDVq1f59re++VHtyjn+qiMEJAHwLdEgAQUIWkoUBBEIwoGwBOEJBFRoCfBd1BgEQkgIguDB24BEIJRD4BHCEmQA4QjB4r0DERBS4L1DKY0QuuU4D5CAUxgTMB6SjsZ4h5IO8IQA0kcokbSEKhxeNICAoNv3cZ8oz0jZLQnWLbdLzknyw4XS8e0ftv3DIsgd4MpDP19ebruPEMI/A/4ZwDPPfO68IfwcfwnOIj710LazaFG8J9rSiBAQISCDb38fFG0UJ+4/07oGHyp05JHKE6gJWIJrWMym5PmMuinwzlI3JUdHhxA8G+ubDPojlJBEKmH14iP4YFGRQOmEqnYkyQgpuoTQRphCtq/rnSGIgJDR/ffVPkDQsu17TwPFOT5afFgE+SzwuBDiBi0x/u+B/8OH9Frn+InAe9LlwrcB3IMN7c3L5U8BgkEEQRulsSTIAAqiCIIw+DDh6PQWR0fvcHp8yOzklOl0TJHPcbaG4JDL9bO1ll6vjwiSuqyRKsER87M//3P0Bl22Lj+Clh1w4IMihE77+gKE9Ahl2qgxLE+7+2Go/cG3Gx5+v+E8gPyI8KEQZAjBCiH+L8Af0V4G/4cQwssfxmud4ycFD5aZgXb52hLHMooMEhGWRHjGJkK32cOwzAWK5fN8RV3uUtV7HB6/xsuvfoPpdJeL66ts9vpsdQJaj1BaEbzHWkcIgoAkijO8l3gfEDLmjZu3efH5P2dWTHn8yadYWb3ChQsfZ239Uygdgdd4B8FXyGiZIngXs4flfj2UM31PxPvg9+f4ceNDy0GGEP4A+IMP6++f4ycRZ3k6seSOcP8m7ldJWgTACYEQHhkaCBYogZLdt77P2298C8KY0Yrk8sBwbTRkEIMwJyAgVilRlOKFwDmJUAlOaIzXeBEThMZ7wSef/hjOW4p6wcnpKbv39tjfPeDjH28YrTxCp7eNimO89eDPco3iXfv+7pt86L0+lJc8x0eCj6xIc45z/K/Bw/Xph5ebgdAWUghoJQneEywEJQgSnK/RosQ1R4wP36Ca3WN68CZpOCTThp7TjGKHAmJbIH3V5g7rGmFLgkwwRDTOI+I+nU4PH3VwIkKICGEheE+/s8pKf4uiXpA3JTv3vsn3X/ozHn30M1zYfor+8Dq4DCECPgSEbG/OWULwaC0J/r3v+hwfNc4J8hx/LdCmDx/OxbURlkC0lWUg0BA8CK0QMlDbHEVJXR+yd+c5Du98h0ScMtANq72aKAQiK5BBo4JA+woVKkIAT0SQASdA6YgkTqm9QjiNjBJ80HinUFYjhUJJ0JEhiRL6YsbpfI+1Vc9rr32Rvd03+fzf+C/p9B6BYJDC45zDG4eONEIojDFoFfF+0qRzfDQ4J8hz/DVBWK4+z0hjmaNbVoCDt4BFRZLga5p6jooK9u6+wvjgDZTZZZRN6DAn9TlxcGgvkTZGuIAIGsmZnEiihGjTlaKtoHsfaExNEAVCRGiVEIQGkRK8wjuPEJJIO2JtyDbWOZntc2EjZTrZ57lv/AEf//TfptPbIkp6KB2jUHhvCF4R6XQZQZ4trf3yvT38ns/x48Y5QZ7jrz4EgFvSRKsLFGcaRw/gkUoSgsPaHKkcSi+49fJXMIs9ovqETM6J3QJl58TCECOQXiAfIqWAwgqNCBKFRISAFh5ramxtkUHjvcGZApl0iDrr1KFLQENw4B3BS7zTBBUxSDaw5Snro5T9o9t87Sv/iqc+/Us8+uRnsU4jRUbwelmUEe+RMZ3pIN/1IZzjx4xzgjzHX3kEWhF428+i2mr1sotFoAnB40yNigPBVxwe3WO8/xp3v/sHXF5LWOnFxK4m9tBNh5hisWxoUYACIQlS4oiwQiBDwAeBQqCDQ1ER+YBWMTY01NUEbIyUDlSfICRSBRQCQkwIgogu48kBK1mXvJkz7C9YlCe88MIXCbHg0WufRAiNVDEETV05kvih01E8XN1evt9z/NhxTpDn+AGE9yzpxHt/+oEV3w9ZCr73fH7vc8TDG94v7yYeUr+0hHE/cmz7Upatfx7vGqSvOT2+zSsvfZPJ7ot8/nJGl5JQL4iDIAoSXweUz0AIvJB4AUFIvPRYIfEiBu/QwRNCQASDJBApCKIhEgolPdaBKwSyfxkvFAKHICCFRoQOvhT0owwhGtJBStYpOTw9pMHy7LNfYjFf8NQTP003y3DOkiSdtji01DyKs86h+wR5jo8C5wR5jnchLOvFZzIaER7qbvYglX7QxxzabWKZKzur0EoFdVUSZylnz7bGtbpAQErwocF7j1IRznmkVIBHCHDOoqTC+YCSmuBASgvSE7zHL0kkSAjeImSBFhNOd17i+NZ3yRa3uLju6LicZBmFKUTbL+0FCAlC3NeOh7OWwPBgsS0JeByNCAgJnva9ITwi0gilUUqQUhC8Q4iaECxVaYAOWbrGfF7S6ytEWdKRhlQHsp5FsuB73/59fF3ymc/8Ctb1SdIYHyQSQfBtq6HwIKRuifM8gPxIcE6Q5/gBhIeiuyAe8lYIHoHAWouSAiklIbQMKhQ462lMQ9bL0FEE3mNsQxQnS4liS0qNMSgdEFISwpK0pMAYRxzLpQQGpJAIIVuJeJDglnpHKWl8wHpHrCuUn5Ifv8Zs59t06jsM4zGDCKRzSH9GeaElQvmwtlDc/6aCeFCkES0hWhHweIKSoDTWh6UhhUQJjRQeij2SKMWHihAsVBWOFN1xDLpgygXaTYgZMzIlTln0eg+bVzz/7S+Ch0999lcxXhNIUTImeE8kVduOLfRy78/xUeCcIM/xbgTxIMf3UNjSlkZ8K7XxBqQDqZDaEYIHJZHOkHUjnJkhVUxVlaRJulwM+zbSkhIlJN56tIqwJhBFGgLEOsaaCh0JrK1QKqZpLFIkaKUgaLyQSOGRWCJKtJtxcPdFdt/8BgNxTC80JDqgxXtUQbTRsVhGjmdbYJnug6XQfPkecaizJW7jEMqjPQQfwFugwYuCRp4gsx6mDkiREaPxqsa7fRqXY+wc7wqEa1A+kGlFXs54ZHuIFhUvfOdPEDrip575dSDDBhCqlf5gzblPxUeMc4I8x7vQxletY454qG1P3F9PB3TkEaKmaUqqYoJ3NVJLqqImSrp0OwOkyMi62bKyWxNpQcCA90vdYIpvHKFpe6ODC6ADSngIHiE8Yrk8lioC2ZKjtQYZShLtcPUx4/1X2L35LGZ2j2xFEXlPLDTONj9gdnrWqPjDsp+CQPAOT0DLNppUUoDz2KrEGosKso1qadsFvQrIoca5CuUTtJQY70hTRVmeULsFUhq0AmVaSZG2sJ5Jps2MJ6+uUpS7PPuNL/D4kx+n20shaIKToEQblp+T40eKc4I8xw9
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAD8CAYAAAAVOD3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9WbBk2XWeCX5r731Gn+4cU0ZOyAQSBECAJEiqOElkdVFUVaukKlnJVP3SD22mp35vPfdTvfVDm3Vb0dq6S9ZtUrUepJJUlDhoAEkQM5ADch5ijht3vj6eaU/9cDwiIxMAKQlIJNi4v1nEdT/u1/1c93P+s/Za//qXxBi5wAUucIELfC/Ux70DF7jABS7wk4oLgrzABS5wgR+AC4K8wAUucIEfgAuCvMAFLnCBH4ALgrzABS5wgR+AC4K8wAUucIEfgI+MIEXkd0TkLRF5V0T+wUf1Phe4wAUu8FFBPgodpIho4G3gvwDuAd8E/vsY4+s/8je7wAUucIGPCB9VBPlLwLsxxhsxxg74n4G/9RG91wUucIELfCQwH9HrXgPuPnb/HvDLP+jJOzs78emnn/6IduUC//+O77cGku95IIJEIBKJ/V1REOXhs/tnyWPP/56fsv5//Ttx/bKP/U4kgETk0W8+fP33YxH5vnv8PXt/gR8jvv3tb5/EGHc/vP2jIsi/ECLy94G/D/Dkk0/yzW9+8+PalQv8JcAHKGVNdA+3x4ckJwBCjKAiqLB+qkRQEMURJBCCJ3gPCAaNipoA+BhRaQJ4YrSABQmAJ8RIjJCoHEiQqBAUrnOYRPd7oiDigQACIUKIQvARJQlGJ2tyDHyYeN+H4n1SvcCPC0qp299v+0dFkPeB64/df2K97RFijL8L/C7AF7/4xYuG8Av8xZDHb/R3JPa3wsMH48NH1jckgiiiQJSECKzqiizLEAIOi6gWTUQFh69aTs+OOJ8e0zRLIg4h0nUdXWsZDCaUxZgiHzLemDDcHOHFIJKgKIhBCF6jxKCUQgNRrSNW7/twU9SHQtyH8aY8dv8CPwn4qAjym8DzIvIMPTH+PeB/9xG91wV+KvB+1NhDPXZX3ueb+HALoCxxTaYRIfQ8RZHnKBXxocX5BcvFAQcP3mN6dJt2fsxyOaNtV2gTKcqMLEtp247ZdA5R4X2/NC/HE4qtLUaTXfZ2nuLa5U+R5zsYPQGniFZABGUAFYmxg6h7goyP9vLR3/BBXvxwzHBBmh8HPhKCjDE6Efk/An8AaOD/GWN87aN4rwv8tOBxgpTHVqjycMsHnyrgAZH1UjxGYgjEECBYquqE5eI+hwfvcOvmd2mbEy5taHaHgcvjiFIaURHEIuKIMRJ2M0xaQNA463ESuHn0Fvvn77F/47vcGL7Mxvgprj/xGa5c+xQmLUFpEFlnPgWRDy+f4/fs9/c+dkGOHxc+shxkjPFfAf/qo3r9C/w0IT72Tx6rpHyIHD/ENWF9eCsCQofEBvyKu7ff5MZb36ZbHTEZCJezFcU4ZWMIufFkqcEkah2kBrwPfV5TKfrrvcE7RWUdu5+6Tt0Ji4Vldn7Ewa0HzI5uUc3usXftecY71zBpgQ9CVAr1fonnA3/D+zsev88fdEGQHxc+tiLNBS7wH4qH5ZjAB7KPECGGSAwRpddL1PCwtKwIUYghoFWHbQ65f/MVTu6/idhTBu6UzbRhJIpBqRkVBSYx+JAQm4BYhSjwwUOMmCRFRCNoovSFlNwEgg9MUuHSTkk9jix2LK2v2b/3p7zx7ld56pM/x/Of+WXyfA9PTiSSED/Ig0AIoNT7edQP5ii/J7S8wI8JFwR5gb90eJwulEBUgu0cWkeEACoiokmjprNLzs5uc/fG11mcvMn2oGGYVxi9wLiWLMIg5iQ+EkOJJwVAglpzcCSgUJIgooiiepKMAAGJNaI9UVtyHShHgQaLaEtuU/bvfYv903v83C/+NuPJdTQpRI3CoJTB2oAgJMmHJckXdcufBFwQ5AX+UiBKX+mNH1pWRwFiJEmk1/ZER/ANEiG0LScP3uP+zReJ7X2e2g4MdIOxcxIqNJ7EK3RQSKMQDWIirOU4IfQiR60FoiPGnhyRgIj0ysZg1vviicqitCfRkO5mFDVkmXDn6C5f+7N/wWd+9td46soLmGREiIEYwzpqTPqI8sNB4vtiygt8TLggyAv8pUCkL3TAo5p1fzs+FGd7vK0wJqJNoFkccvPlr3B68B7D3LMzAW2nlNqRGkGRIGhEFBI0RE8MFUQLktDrESMiBkTjQwRUX2h5VHgxiMoADUEjUaHX8W0InjIqzNAQHTw4O+TN7/wB9TMnPP/JX6YsJr3uUqcoEWznSRLDBSv+ZOGCIC/wlwgfLm/EvqgRHTF2KDqiRM5PH/Dmi19iefvLfPrpXQwe3XVsbY3w9QJNjuisF4Er8H2akRA6cB2aBFEGEQ0xEoMnhD5v2C+xVf+YgNIZEdVHkkEgBpQDHQPjJMWKQAGDvZS7B4e8/epXCOR85jNfJNFFL/nBoM33CR8vePJjxwVBXuA/Hf/BRdYfdKb/hxYeBInrzOO6ePEoDylADIjuW2fOHtzk9de/w+m9V/nMnmekDhgNJkjQxMWcBINICgFC9FgcXgJBIsobxCq805jEoHVCCLoXoSu9JkgB0UjURJ1gJfaEqTJ0EGTdvZMaTdc0DIqENIPUzZBJxJ5XvPP2izhnee65zzAeX0LEICrlcd+Y9xMJFyz5ceKCIC/wQTwmvn6sLeX9h0O/3H1UcI0QvEdp+pwdoJV+9JgoaOqWLDfrHFsghIiSfjnpQ0BJ3x4Yoc/tKcVDPgyhzzqqR8QYsN6DMijREBVaKQgNs/3XuPfmlzGzO3z2smcYOnISaCqCjcQgKAyC6nWNRMATcTjvMSElIYWgwBuiKKLqyVNUJEQQrVCJEKLvu29cgohBrTsE3fpVVVxrL9uWNIlspwpVd9Sl54G9x803DrHVKZ/5/F9lOLpOlL4BUa3/SYgoEQh+raX8iL7vC/y5uCDIC/w5+OAyz/u+qKDWxBU8KMUj4hMBEU0MvdomhoDSgtEaQXBdg+ARpQjSS1tE+ugseLdebvbV3BgjIr2wWgB8L98RExDxdCGiRJGIEJxltv8md179Q3T9Hk+VHapdMSyGRAveWZKHnTcxELoOYkS0ItWCIaCCR5OQqAS0Bh0JsSM4S1ABrRQheDQGbTJ81yEqR7sB6AiJEHHY6PBBSJMSnQo4R6gbxNdsJoqQWGg7xAr3b36Tohzyuc/vEkxGkHWOEyCsL1Dh0UdygY8BFwR5gR+A79XePWwjhtAHNSoAgRAddV2RpgVGp/gIMSi0VigFMRi8E5QqELWW4ogmBLtm14BIROk+32etJUZBKY1SCusCGtVHmji0ipgQ0AFUaDi88wrvfPdfspnc58mrOWmXYNIB3lW9SDxKr2f0HokR7zwPCzCiDdpocqN7wqYDo0EiyjX4roZgUcZgnEOiIEmCai1iCoIKqCwHL0SzlqYrDcERUXjvcW2LDpa8TBjolO1EgUTsUcNb3/0KkPOZL/wW2ow/0IIYQ0T0xSn6ceLi07/AByEfakd5jCSVFry3oDoES9MumS7Oe1cbGzAmRTCMR1vk2RCIxKjwPmCM7iPHEHHOkaU94REDrnMkaUIIHhEhekuaF3jnUUqTGkNUQogR5wVFJNMG3yw4vP1d3nvzD1HtTSbjDhMSfKNQwaAlALbPN4ZAdJ7oH7YRqp6bRRAla9syB9H30WrwYFt01yBdR3Ae5SMqrq3LqhZnCuIwYMoSlSkk1yRZjlYKa1tCNCgRkixDPLS+I8aEvdEW1fKEyxsJ90+nvPKdP+D69SfZ2vs0IZSEAEYJPkaMulhbf5y4IMgLfAgP7biEvq3ufYQQiHRo6ZgtD3jjrRc5nx6RJAUxaIgG7+FnP/vzbG9dIU0GtJ3n+PiMre1tsqRAlMIYQ9s0RDx5nuGsXecZ+/xd267I8hRnO5QySASPwgtEyVAxIK7h9P6L3Hz99ynlHpf3FDkebTVpmlEvlqS6QXDrPyn0+dO4ThWg0VqvI1iBGFH4XkfZNoS2IXYdOIt0jlB
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "path_image = \"images/fruits.jpg\"\n",
    "#for i in range(1,13):\n",
    "    # Pc, index, img_seg = kmeans_image(path_image=path_image, K=i)\n",
    "plt.figure()\n",
    "plt.imshow(io.imread(\"images/fruits_2.jpg\"))\n",
    "plt.figure()\n",
    "plt.imshow(io.imread(\"images/fruits_4.jpg\"))\n",
    "plt.figure()\n",
    "plt.imshow(io.imread(\"images/fruits_10.jpg\"))\n",
    "plt.figure()\n",
    "plt.imshow(io.imread(\"images/fruits_255_cuda.jpg\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La clusterisation avec 255 couleurs à été réalisée à l'aide du script `Kmeans_skcuda.py` afin d'accélérer le temps d'exécution. \n",
    "\n",
    "Nous avons vu en A4 la puissance de CUDA pour le traitement de données importantes, ce qui est notre cas avec le nombre de pixels de l'image et le nombre d'itérations qui augmentent en fonction du nombre de clusters recherchés. J'ai ainsi décidé - pour voir à partir de combien de niveaux de couleurs peut-on apercevoir une image nette - d'observer le traitement Kmeans jusqu'à 255 clusters.\n",
    "\n",
    "On observe ainsi qu'au dessus de 40 clusters on arrive à bien distinguer les fruits et légumes présents sur l'image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEWCAYAAAAtuzN2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABHQklEQVR4nO2dd5hU5fXHP4elL0gHadJVsACKvRGxoSJqFLEFS+wa/UWjxjRjYgzGrtEEu7EgFkSNDRBUFEEB6SK9L2XpnV3e3x/n3ty7y87s7M7MTjuf55nnvW3ufe+075zznvcccc5hGIZhGOlCtVR3wDAMwzDCmDAZhmEYaYUJk2EYhpFWmDAZhmEYaYUJk2EYhpFWmDAZhmEYaYUJk5EyRKS3iCxL0rmricgIEbk2See/R0SeS8a5swkR2U9EtohIXqr7YmQOJkwZjvel9x97RGR7aP3SVPevIojIIhFpn6DT/RX43Dk3JN4TlSWgzrm/Oed+Ge+5y7jWFSIyLrS+j4h8LSLviEjNRF8v2Tjnljjn6jnnilPdFyNzqJ7qDhjx4Zyr5y+LyCLgl865UanrUXrgnLsn1X2IFxFpBHwGzAV+4ZwrSnGXDKNKMIspS/FcWXeLyHwRKRSRYSLS2NvXXkSciFwpIktFZL2IXC8iR4jINBHZICJPhc51hfev/SkR2SgiP4pIn1L7F4jIZhFZGMlSE5E6IvKSd71ZwBFR+l9LRB4SkSUiskpE/iUidbx9H4nIw6Fjh4rIC95yAxF5XkRWishyEflr2I0kIteIyGyvr7NE5DBvuxORzqHjXvKemw98DLQKWaKtROReEXk1dPw5IjLTe+3GikjX0L5FInKH99puFJE3RaR2Oe9fM2AMMAO4zBclr19Pi8jHXl++FpF9ReQx73X9UUR6hs7TyrO21njvza9C+44UkfFen1d672/N0H7nfS7mesf8U0TE29dZRL7w7metiLwZ4T78z1p1b32s97p+4/X/AxFpIiKvicgmEflOQlaziDzufUY3icgkETkhtK+OiLzs3fdsEblTQpZttHs30hznnD2y5AEsAk7xlm8FvgXaALWAfwNvePvaAw74F1AbOA3YAbwHNAdaA6uBk7zjrwCKgP8DagAXARuBxkA+sAk4wDu2JXBQhP79HfjKe15b9Ed3WYRjHwXe946tD3wAPODt29fr38nApcACoL63b7h3r/nevUwErvP2XQgsRwVRgM5AO2+fAzqHrv8S8FdvuXfpfgL3Aq96y/sDW4FTvdfnTmAeUDP0vkwEWnn3Mxu4PsJ9XwHMAmYCzwBSav9LwFrgcO+9+xxYCPwCyENdmGO8Y6sBk4A/AjWBjt5rdbq3/3DgaNRz0t7r122hazngQ6AhsB+wBjjD2/cG8DvvGrWB4yPcT3vvPNW99bHea9MJaODd60/AKV4/XgFeDD3/MqCJt+92oACoHfo8fQE0Qj/n0/z3qbx7t0d6P1LeAXsk8M0sKUyzgT6hfS2B3aEfIQe0Du0vBC4Krb/j/0h5P5Yrwj+S3g/t5agAbAB+DtQpp38L/B82b/1ayhAmVDS2Ap1C244BFobWfw4s9X6kj/e2tQB2hvsBXBz6of4UuDVC3+IRpj8Aw0L7qqEC2Dv0vlwW2v8g8K8I/bgC2Oy9V0eVsf8l4NnQ+i3A7ND6IcAGb/koYEmp5/+W0A9/qX23AcNLvSbHh9aHAXd7y68AQ4A25bzn/mctLEy/C+1/GPg4tN4P+CHK+dYD3UOfp9ND+35JIEwVund7pNfDXHnZSztguOeC2YAKVTH64+2zKrS8vYz1eqH15c77dnssBlo557aiFtT1wEoR+a+IHBihT61QMQmfoyyaAXWBSaH+f+Jt9/kAtRDmOOf8YIF2qMWyMvS8f6OWE6iVNj/CNeOhFaF7cc7tQe+zdeiYgtDyNkq+tqWZCtwBfBx2y4WI9X1rh7ogN4Rej3vwPgMisr+IfCgiBSKyCfgb0LTUtSL1+070D8REz4V5VZT7qWz/8Vygsz2X4QbUyvL7WPrzFF6Oeu9GemPClL0sBfo65xqGHrWdc8sreb7W/viCx36oFYVz7lPn3KmoVfYj8GyEc6xExSF8jrJYi/5AHRTqewMXCvQA7kfFtqWIXOxtW4paTE1Dz9vHOXdQaH+nCNfchoqhz76h5fJS8K9AfwgB8F6ntqjVVCmcc4+jrqqRInJwJU+zFLUyw5+B+s65M739z6DvVxfn3D7oD7dEOlmp/hU4565xzrUCrgOeDo/RJQJvPOlOYADQyDnXEHUh+31cibrwfMKfrfLu3UhjTJiyl38B94tIO9DBdBHpH8f5mgO/EpEaInIh0BX4SERaiEh/L0hgJ7AF2BPhHMOA34pIIxFpg7qh9sKzOJ4FHhWR5l7/W4vI6d7yicCV6LjKIOBJEWntnFuJRrE9LBpmXU1EOonISd6pnwPuEJHDRensvz7AD8AlIpInImcA/nNA/9E3EZEGUe7rLBHpIyI10LGQncA3EY6PCefcg8DjwCgROaASp5gIbBaRu7xAgTwROVhE/KCT+uj44BbPyr0h1hOLyIXeewjqXnNEft8rS310bHMNUF1E/gjsE9of/jy1Bm4O7Svv3o00xoQpe3kcDR74TEQ2o4EQR8VxvglAF9SauR+4wDlXiH6Gfo1aDevQH/RIP3B/Rl1eC1EB+U+U692FDpJ/67mZRgEHiMg+6PjGzc655c65r4DngRc9S+UX6GD3LPQH823UksM595bX99fRcZz30GAE0GCRfuh42aXePrzn/YgO9i/w3EKtwh11zs1BB+mf9F6ffkA/59yuKPcXE865v6CCOlpEIll7kZ5bDJwN9EBf87XeuXyBvQO4BH0tngXKjKyLwBHABBHZgn7ObnXOLahI/2LgU9SF+xP6udlBSXfdfcAy9N5Goe/1Tojp3o00RkoOGxjG3ojIFej8qONT3RfDiISI3AAMdM6dVO7BRlpjFpNhGBmJiLQUkeM8l+0BqAt1eKr7ZcSPZX4wDCNTqYlGXXZAXbBDgadT2SEjMZgrzzAMw0grzJVnGIZhpBUZ7cqrVq2aq1OnTqq7YRiGkVFs27bNOefS1jDJaGGqU6cOW7duTXU3DMMwMgoR2Z7qPkQjbRXTMAzDyE1MmAzDMIy0woTJMAzDSCtMmAzDMIySiByAyA+hxyZEbkOkMSIjEZnrtY2ScvlMnseUn5/vLPjBMAyjYojINudcfowH56GZ8o8CbgLW4dzfEbkbaIRzdyW6fxlpMYlIPxEZUlxcnOquGIZhZDt9gPk4txjoD7zsbX8ZODcZFzSLyTAMI8cQkV3A9NCmIc65IREOfgGYjHNPIbIBrYsFms1//f/WE0hGWkypwjl48UXYFXcxA8MwjEpSWAh33w1z5sRzliLnXK/QI5Io1QTOAd7aa59aNUmxbEyYKsCUKXDVVTByZKp7YhhGxjNlCrz+esWfN2ECDB4MBQXlHxs/fVFraZW3vgqRlgBeuzoZFzVhqgDbvbnS5j00DCNuHnsMbimjiPP8+TBmTMlt//0vfPmlLn/3HYjAYYclvYvAxWiRTJ/30arReO2IZFw0o1MSVTU7d2q7Y0dq+2EYRopYtQqqVYNmzeI/19q1sGED7Nmj5/R54AEYOhTWr4caNXTb2WdrW1wMEydC165Qv378fYiGSD5wKnBdaOvfgWGIXI1WFR6QjEubxVQB/LElEybDyFEuvxyuvTYx51q7VkVpy5aS2wsL1S0zderezxk5Ui2mI49MTB+i4dxWnGuCcxtD2wpxrg/OdcG5U3BuXTIubcJUAUyYDCPHWboUVq5MzLnWrtV2w4aS2/31ceOCbXXravv738OaNXDEEYnpQ5piwlQBfGHyXXqGYeQYGzbA5s2JOZcvTOvX730NCITJOf3xycuD77/XbVkuTBk5xiQi/YB+tWrVqtLrmsVkGDnOhg3BuE887NoFmzYF5yx9DVBhck7/CRcVwd/+Bp07w+efQ48e8fchjclIi8k594Fz7tq8vLwquBbceitccEFswlRYqG5jwzAyhMmT4Z//LP+4HTv0UXpMqDIUFgbLZQlT8+Y
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cpu = []\n",
    "cuda = []\n",
    "i = 0\n",
    "with open(\"timing.txt\",'r') as data_file:\n",
    "    for line in data_file:\n",
    "        data = line.split()\n",
    "        if(len(data) < 4):\n",
    "            i=i+1\n",
    "        if(len(data) > 3):\n",
    "            if(i == 1):\n",
    "                cpu.append(float(data[2]))\n",
    "            if(i == 2):\n",
    "                cuda.append(float(data[2]))\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_xlabel(\"K\")\n",
    "ax1.set_ylabel(\"temps (s)\")\n",
    "ax1.set_yscale('log')\n",
    "ax1.plot(cpu,'b')\n",
    "ax1.tick_params(axis ='y', labelcolor = 'blue') \n",
    "plt.legend(['CPU'],loc='lower left')\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.set_ylabel(\"temps (s)\")\n",
    "ax2.plot(cuda,'r')\n",
    "ax2.tick_params(axis ='y', labelcolor = 'red') \n",
    "\n",
    "plt.legend(['CUDA'],loc='lower right')\n",
    "plt.title(\"Temps d'exécution Kmeans image\")\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.4 64-bit",
   "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.9.4"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "2ef431f6525756fa8a44688585fa332ef3b2e5fcfe8fe75df35bbf7028a8b511"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}