{
 "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": 188,
   "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": 189,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean = [1,2,3,4]\n",
    "# sd = [0.25, 0.25, 0.1, 0.2]\n",
    "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)\n",
    "path_image = \"fruits.jpg\"\n",
    "# print(mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fonctions à utiliser pour le clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_points(mean=1,sd=0.5, nb=100, dim=2, clusters=2):\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": 191,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance(points,Pc):   \n",
    "    return scipy.spatial.distance.cdist(points[:,:], Pc[:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kmeans(points = [0,0], K = 1):\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": 193,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors=['red', 'green','yellow','blue','purple', 'orange']\n",
    "def visualisation(points, index, Pc=[0,0], K=1):\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],'r+')\n",
    "        plt.grid(True)\n",
    "        plt.axis([min(mean)-1,max(mean)+1,min(mean)-1,max(mean)+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [],
   "source": [
    "def img_2_mat(my_img):\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": 195,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mat_2_img(mat,my_img):\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": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kmeans_image(path_image, K):\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [],
   "source": [
    "points, mean = gen_points(mean,sd,nb,dim,clusters)\n",
    "# print(points.shape)\n",
    "# print(points.mean(axis=0))\n",
    "# print(points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist = distance(points,points)\n",
    "# print(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pc, index, clusters = kmeans(points,K=K)\n",
    "# print(index)\n",
    "# print(clusters)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtwklEQVR4nO2de2wd133nvz/SJEWaDg3RIteIw8t4Ywk1wjSFjawWAVbPArZibxGjBWRcCUpSgRurLqSFhWwDAhUF9GaxhVBIwMYOVKWwIBEVBNRtaltFEOoVLFBtY7WOadeRYyikNkUhQhQimxZNSeTZP+YOOZx7zjzPmZkz9/cBLsg7d2bOGT6+85vf65AQAgzDMIy9tOQ9AYZhGCYdLOQMwzCWw0LOMAxjOSzkDMMwlsNCzjAMYzn35THogw8+KL7whS/kMXSmfPLJJ7j//vvznoZxmuU6gea5Vr7OYnL58uUbQog1/u25CHl/fz/eeuutPIbOlAsXLmDjxo15T8M4zXKdQPNcK19nMSGiKdl2dq0wDJOasYkxDB4eRMvBFgweHsTYxFjeU2oqWMgZhknF2MQYhl8fxtStKQgITN2awvDrw7mKebPdWFjIGYZJxcjZEdy+e3vFttt3b2Pk7Egu8ynijcU0LOQMw6Ti2q1rsbabpmg3lixgIWcYJhUDPQOh27N0dRTtxpIFLOQMw6SitqWGrrauFdu62rpQ21IDkL2rI8qNpWywkDMMk4rqUBVHnz2KSk8FBEKlp4Kjzx5FdagKIHtXR9iNpYzkkkfOMEy5qA5Vl4TbT9auDu8N5NqtaxjoGUBtS005vzLAQs4wjFEGegYwdauxjsWkqyPoxlJG2LXCMIxRmtHVkTUs5AzDGCXMh86kh10rDMMYJ6mrY2xirKl83UlhIWcYppC4aYtuxoubtgiAxdwHu1YYhikkzVihmRQWcoZhCkkzVmgmhYWcYZhC0owVmklhIWcYppDknbZoUytcFnKGYRJhWujyTFu0rRUuZ60wDBObrDJK0lRopkldDAq0FjFjhi1yhmFiU/SMkj1v7sHO13YmtqhtC7SykDMME5siC93YxBh+8NYPICBWbI9zo7Et0MpCzjBMIDJfuA6h0+Vj959n7z/sbRBxl6g3mrwDrXFhIWcYRokq6LftsW2phE5XMFF2npm5GeX+AiLSTcO2/jDahJyIWonoX4joDV3nZBgmX1S+8DO/PJNK6HT52GXnCSPqTaM6VMXkvkksHljE5L5J7SKuM+tHZ9bKXgDvA/iMxnMyDJMjQb7wNBklunzsSX3y7k3j1S+/muj4tOjO+tFikRPRIwC+BuCYjvMxDFMMTAX9dJ1XtX9vZ+/S04KKPAOzurN+dLlWDgP4DoBFTedjGKYAmAr6JTmvzBWhOs+Rp48suUUqPRXp+fLMQNGd9UNCyKO7kU9A9AyAbUKIPUS0EcB+IcQzkv2GAQwDwJo1a544ffp0qnFtYHZ2Ft3d3XlPwzjNcp1A81yr9zrHr4/j2K+OYXp+Gn0dfdj9+d3Y2r819Rj+865fvR6Xbl6SjjN+fRyHPjiE+cX5peM7Wjqwf+1+AAicX9Cx6+9fn8vvc/ul7bg+f71he39HP06tP6U8btOmTZeFEE/6t+sQ8v8JYCeAewBWwfGRvyaE2KE6Zt26deLKlSupxrWBCxcuYOPGjXlPwzjNcp1A81xr1tfp9xkDjmXtBlAHDw9K1/2s9FQwuW8y0vllVZ55/T7DrlcFEUmFPHWwUwjxXQDfrQ+yEY5FrhRxhmEYP2El8XFdEXve3IOjl49iQSyglVox/MRwJMHPClesda1+xL1WGIbJHZUgu1b4QM+A1CKX+bn3vLkHr7z1ytL7BbGAV956BR/MfIAPb36Ia7euYXXnagDAzbmbGHg7nyXk0mT9+NFaECSEuCDzjzMMU1505EOrAo8ECgxqyoKjRy8flZ7r7K/OLhUOzczNYGZuJnIxUtFb2nJlJ8MwiZFVVu54bQce+vOHpGKnEsTalpo0VVBALLlXohYgLYiF2NcRlPpnQ0tbFnKGaWLSWpqqysqZuZkGsQsSxOpQNbQ/StRKy1ZqjXUN/nH8FL3TI8BCzjQJE2MTODx4GAdbDuLw4GFMjE3kPaXc0WFpBuU9+8UuTBBV+d5R+6O4DD8xHGk/P67f3E+ROz26sJAzpWdibAKvD7+OW1O3AAHcmrqF14dfjy3mE2MTuLT9UmluBjoszbCiGq/YhQmizA/uEucm8/LXXsYLT76wZJkntdBdbGhpy0LOlJ6zI2dx9/bdFdvu3r6LsyNnI5/DvRnMX59PdTMoEjoszSDxBVaKXZggev3gMuLcZF7+2su496f3IA4I3PvTe8pzerk5d1O63YaWtizkTOm5de1WrO0ydNwMioYOS9MV397O3obP/GIXJojeoh0Vys9GRwPnGXbDAdTXbUNLWxZypvT0DPTE2i5Dx82gaOiyNKtDVdz4zg2cfO5koNgFCaLfX69CeZM5eDB0jkHWPgBse2xb4PEmW9qmhYWcKT1balvQ1tW2YltbVxu21LZEPoeOm0HR0G1pphG7KH3Fo95kVJk47vxUYn7ml2ciz7docGUnU1omxiZwduQsbl27hc7Vnbiv8z7M3ZxDz0APttS2YKg6FPlcW2pb8Prw6yvcK3FvBkVEZ3VhGEE9uIPcKQSSl7CPjq60xMnJQ5/cfB+m/su9hjHcY23IQokLCzlTStzgpCu8czNzaOtqw3Mnnosl4C7uMW++9Cbmp+cT3QzKgqoBVRhBWTKqEvzAplijo8u+cSJACGlzLW+QdOTsiNJ100ItaDnYsnRN7v46eqGYhoWcKSVBwcmk4jtUHcLMZ2eaovuhijQr2wRZwieeO4Fv/t03cXfR88TT0tbgSgm7iQT1bPF3G/TjVoRO3ZrCt370LQghluaTdgUf07CPnCklZQxOFoE0uedhWTJEK0v0/e8DC5gOHAgco5ValSIuyzO/s3BnxU0FKF41pxcWcqaUlDE4WQTS+JeDsmRGzo7gzsKdFZ/dWbgTvTK07mJRjaHqv0KgWL1ZiupHZyFnSomOTJWyoLNzX5rc86AsmSg3iCj7qMZQZar0dfTFqvwsUjWnF/aRM6XE9YO7WSvNGpzUvVp7bUtNurJN1NxzVZaMKtjp7X8StSe5agzZvHd/fjdqv4g296JVc3phi5wpLUPVIeyb3IcDiwewb3IfhqpDTdc8S3fnPlNVjrUtNbS1tDVs//jOxyta3SYtYFLNe2v/VqW13t3eHfs68+pbzhY50zT4UxLdfikASmupm8iZNpF7Xh2qYu8/7MXM3MyK7a6f3Dtm0pRA2bwvXLiA2pYavvWjbzX46Ofvzcc6v+6nnziwRc40DWXslxKGDZ37XFRNq/w+cLd61A2SprV+q0NVPND+QMP2u4t3Yz25qJ5+vv3Gt41b6SzkTNPQjCmJNnTuc1HdXFqopUH8dK3aM359HIOHBxueBFziPLmo9p29M2t8dSEWcqZpaNaUxM77Ope+7+3sTe3TNuUHVnUoXBALK8RvbGIMu/52V2rf/9jEGA59cEgaQHWJ8+QSdV8T+egs5EzT0Gwpia7V6rU25+7NaTmnCQvTDUjK0gFv372NHa/tQPf3uvHNv/umMvc7jgU9cnYE84vzys/jPrnE2Vd3PjoLOdM0DFWH8OzRZ9FT6QEI6Kn04Nmjz5Y20GlirUnT61dWh6pYFIvKzz+5+0lDxaWXOBZ0mJh6n2SiUB2qSvuyy9Ado+CsFaapGKoOlVa4/ZjIWMmic6AqXzyMKBa0t1dLC7UEVnW6C0gD0bNOjjx9JLSnS3tru/YYBVvkDFNSTGSsZJEFE2U1Hz+t1Brq+/e7haKU5sd92oiygMUD7Q9oT0dkIWeYkmIiYyWLLJggX7mMrrYuHP/68VBxVC1e0UqtIJDkCIe4TwduiqTqnKo0yzSwkDNMSTFRhZnV+pXVoSqOf/042lvbGz5rpVb0dvbGHl/l/lkUi1g8sKi0ogmUKJibZQ4/+8gZpsSYqsLMoie3O4a34rO3sxdHnj6SaPywXi21LTXsfG1nw8ITAmKpujQOafvSxIEtcoZhCou7sLM4ICAOCNz4zo3ENxGZW6ijpWNJWKtDVeXqQUmCuVk9vQBskTMMU0CSLicXhKxXy46Hd6w4b6WnEqnDYpwxs3h6SW2RE9EqIvonIvo5Eb1HRAfDj2IYpsi4petZd/EDzBcdub1aJvdNYmv/1hWf29TSwIsO18o8gM1CiN8G8GUATxHReg3nZRglzdaONku8pesm+4OoMF10FESW7hCdpBZy4TBbf9tWf8kdTQyjAbcd7a2pW4BYbkdbZDEfGwMGB4GWFufrWHYGbmxkpetZrleZRdFREH6rvegiDgAkRHrNJaJWAJcBfAHA94UQ/0OyzzCAYQBYs2bNE6dPn049btGZnZ1Fd3d33tMwTtbXeWn7Jcxfb+yR0dHfgfWnzD4MJrnW8fE+HDq0DvPzy3nRHR0L2L//CrZunQ499tixRzE93YG+vnns3n019Ji0bL64WRr0IxDObThndGwA2H5pO67PX2/Y3t/Rj1PrT2kdy7b/0U2bNl0WQjzp365FyJdORvQggL8F8MdCiHdV+61bt05cuXJF27hF5cKFC9i4cWPe0zBO1td5sOWg/JmPgAOLB4yOneRaBweBKUlNSWsrsLgIDAwAtRpQ9Rl+Y2PA8DBw2+Nl6OoCjh5t3Fcng4cHpQG/Sk8Fk/smzQ1cx79AA+D4qU24OGz7HyUiqZBrTT8UQvwGwHkAT+k8L8N4sa0d7TWFR2BhARDCEfnh4UZ3y8jIShEHnPcjdQ+HKXdNbUsNHS0dK7ZlGfCz1U+dJzqyVtbULXEQUSeA3wXwi7TnZZqA0dGlb+MEL21rRzsQIXPNK9AuqhvAtWvL1vrUVPDNIAnVoSr2r92fq5Da6KfOEx0W+cMAzhPROwB+BuAnQog3NJyXKTsHnUzVuMFL29rR1mqOSyQMv3CrbgBCALt2BVvrYYRZ81v7t7KQWkTqgiAhxDsAfkfDXJgQJsYmcHbkLG5du4WegR5sqW0prHjFIWgtTdX12dSO1vVnj4w4Yt3S4rhV/PiFu1Zr9JG7yI4H1Fa8F7/v3bXmvXNl7IJL9C3BxpQ7KaOjAJHzAgAi7Jv679iA8w27lmktzWoVmJx0gpvHjzda6F1djnD7jzl6FOiNtlYBAGD16nC/eZjvXQemloNj5LCQW0JpVoAfHXV8A262lBD4897/hYvY1LBr5+p4K7TYgivQlYpzP6tUgjNR5iKuztbeDnz0UbjfPMj3rgOTlZmmsfUGxEJuCc24AnyZ8Vrok5NqEZdZz15aW5dvBg88ANz1rYIms7RVvvcoQdko5FmZmQabb0As5JZgW8pdJA44Od9zN+Ump2q7DfiDiXv2JEsVDLKSu7ocN417M7ipWK/Afw5Z8FXm2klK3pWZSbH1BgSwkFuDbSl3kainH5btJiVLDXzllWSpgkFWcqfP8xTV0o7r2olLlgsq6MTWGxDAQm4NtqXcxaFsN6kwdwgQPbgYlLo4M7PyhrBt23IM2UVlaUd17STB1g6Ctt6AAO5HbhU2pdzFwb2msqRWRg0aRtnPm7ooK/P33hCOH1+OIQOOqO/alX1Koazvt45+4qbJckUf3bCQM4WgTDepgQG56Mr2i0K16rxaWlYKtcu1a/KnACGAM2eijTE2tpznPjAA7NjRhzQtSLJaUEEntt6AAHatMIx2olRyel0eUXumBPnA06QUynz6hw6tK3SrXVMUsTWANyUS/ZBaOyzkDKMZWTDxhRfkwcU4PVOCsk3SpBTKrPn5+VatBUJMMvwpkWhBu2w/dq0wUsraDiArXHdIGEFVlv7j/aX+/va3spa3UVIKTRcIMcmRpUTKYIu8ZOhYAq007QA0MzYGbN++Xmvb2Lgiqso2qVadwGZrfe2K1tbogU7TBUJMcqKmPrKQlwhdAlyadgAacV0g16+v0to2VpeIjo05WStuM62FBed9lPnJXDYdHQvaCoSY5ERNfWQhLxG6BJjbATRiqtGUrirLNPOT+fT377/CnRALgCwnXwYLeYnQJcBlq7TUgSk/sq4qy7Tz87tsTKwLatMC1EXBv1oSFnFHth8LeYnQJcBlq7TUgUk/so4qy6L7uU2uaFR2vCmRuA6pn5SFvEToEuAytwNIiulGU2kp+vyy6IHezHD6YYnQWerurbR0UxFf2/la06YiulbySy99iunpVQ2pf3kTlpqYN5ziaBYW8pKhu9TdzYRxg6huJow7Vtnxl67v3n0Vf/Znj+c9LSlRc9fzQNW2oCiuH9th1woTSDOnInLpuj6K7vqxHRZyJpBmTkXk0nV9mO6B3uywa4UJpGegxykwkmz3lvG/PfB26Xzn7NfVS5FdP7bDFjkTiCoT5rFtj5W+jL/oKX0M48JCzgSiSkX85Zlflt53zqXrjC2wa4UJRZYJ89rO16T7lsl3Lkvp27HjCqpVs1kr/kyZIqURMsWELXImEc1Sxp9F6bqX2BWQ9QWsmeaGhZxJBJfxmyF2BeTBg8bnxBQfFnImEVzGb4Y8MmW4mZX9sJAziRmqDmHf5D5sOLcB+yb3sYhrIFKmzOiok4xN5Lx3v0/gZlG5csbH+2Kfi8mP1EJORJ8jovNE9K9E9B4R7dUxMYbJmzws1UgVkKOjjuoK4bx3v08g5CpXzrFjj8Y+F5MfOizyewBeEkI8DmA9gD8iomI2o2CYiCS1VNOKf9YVkCqXzfR0h5kBGSOkFnIhxL8LIf65/v3HAN4H8Nm052WYPEliqerquR2rP/mBA/FO7kPlyunrm091XiZbSLiPZzpORjQI4KcAviiE+Mj32TCAYQBYs2bNE6dPn9Y2blGZnZ1Fd3d33tMwThmvc/PmDRCCGrYTCZw7d1F6zPbt63H9+qqG7f39n+LUqUva56iD8fE+HDq0DvPzrUvbOjoW8OKLP8czz3wUcGQ5sO1vd9OmTZeFEE82fCCE0PIC0A3gMoDnwvZdu3atMMtJIURFCEH1rycNjyfn/PnzGs5SjGsJQs91FotKxXU8r3z1988pjyGSH0OkZ04nTzrzInK+ntT0pyA7bxl/pzJsu04AbwmJpmqp7CSiNgB/A2BMCCEv+cuMMTiGv/tcPFV/DwC2lccV81q8zbJ6Bnrw8I6HgY25TccItZrjFvG6V7q6nH7kgDwEZLLntuu2cefjum2A9P5zWTOrCxfSnZPJFh1ZKwTghwDeF0L8RfoppWUEy8Lncru+3TaKdy3uQhPeZlkfHPqgVM2yAHXQMaiy02TPbV4qjQlCR9bKVwHsBLCZiN6uv7ZpOG9CVJUTNvYeLd61yBaaWJxfLFWzLJe4iyKbzDjhlrpMEKldK0KI/wOgMSqUGwNwXBCy7bZRvGtp5oUmomCq5zYvlcYEUcLKzhoA3/MtuurbTTIGYBDOj3Sw/j4teV2LmmZpllU0irpUWp7l/dxaYJkSCnkVwFEAFTgPCpX6e5PBQTcoOQVAwA1K9vWNpzxvHtcSxBi21H6Etq47K7a2dLRwsyzDFHGpNF1587aNXURKKOSAI3STABbrX03/tcuDko8+ekzDubO+FhXOzWqoegHPHn0dPZXfACTQU7kPa/ev5T4rGRDXZ+9iynLNMwDLwd+V8MISWpBHnDo6zPauzpblm9VQdQJDVTdLpReffnofgOfh+O5rsC/Ns7yYTFvMMwDLwd+VlNQizxp5xGl+vkwd5FT/ITNYteo6vC4lPfGBaLCfNBiTlmuea5ryeqorYSHXgjwoefXq7jwmY4io/yHZ5bnr8JOW/UZg0nLNMwBb1OBvXrCQx0aWnSIPSk5Pb81nikaQ3axUZPN8m9babIaAmUnLNc8AbBGDv3nCQh4LeXbKsphPIrugpIl0xyBkN6texb5mn29dK1qWVw1EtzabIWBm2nJNGoC1feyiwUIei6KUzAfdUNKccxDBNwb/zeoIss5z91rRKqJam2ULmMncRGy5NgectRKLopTMB91QkvyHJm3O5Xz26acvYdWqaWSRtSKzor3EsTbLVC0Zlp3Cwl1u2CKPheo/POv/fN03lDRPGlVcunQK2biUxnDhwiAWFlrwq18N4vnnVz41xLU2yxQwM+kmKntAuAywkMdCFvAjOBbsILJLu9N9QynKk0YQzlPD4OAUWloEBgen8Jd/Obwk5pVKfD9pmdwOKndQkAsqCs0QEC4DLOSx8Ab8AEfE3RWWssyh1t2DpShPGkE0PjXcf/9tfO97I6ms6LIEzFTuIKJ0otsMAeEywEIeGzfgV8GyiLtkFfjU3YOleM25GpGbnAMD16y1onVSqzmi7UeIdKJbtoBwGLa6kVjIE5O3O0KW7uhmnhCcODZB7vLxZ6gAxWrOJUNucra0DDS9iAPOjUy1/G4a0U2Th26bKNrsRmIhT0zR3BHelEQAWKh/9bt8VKmLQDGac6lQrVWS4xomBaNSkW9Pk4WTNCBsoyja7EZiIU9M0dwRsswTF6/Lpyi58HE5E3N782EiCydpQNhGUbTZjcRCnpii9QoP+2u75vsa9/i8sXXe5nFdGDt3Ap2dQG+v3iycJAFhG0XR5kZc5RLyzJ1yRekVDoS7dAZ8X+Menze2ztssfhfGzAwwNwecOJFvFo6NomhzXUF5hNxGp5xWgppaeV0+UV1Cpnu5xD1/0VxZxaCoLgwbRdHqugIhROavtWvXCu1UKkI4Er7yVanoH0vJSSFERQhBQoiKeO+9kQzH9o4PIURr/Wulvl22HwV83iVW/tq6JPs5nD9/PsE8o58/+rzNE/9azUIk/7MnSndeHdd58qTz70fkfD2Z/a8rlKL9PsMA8JaQaGp5eq3k7pRr7Feybt0hAL+F7Fwu1Yhjhe2nu5dL1PPv8sxPRtTrax6K3C+Ge7xkR3lcK6q/3JaWjNwrjeLU2jqPeNkgWbemVWE6sKg6zwKyXmHIdmx0YTD6KY+Qy/6iAWBhISNfeVrxM9GaNimmA4tB57EhFbI4WO3XZbRRHiF3/6JbWxs/yyT6E1f8/Nb3XgTnd2dprZsOLNbgpGyqKHCOWgHJol+MbVWazUZ5hBxw/oIXF+WfGfeVN4rfwkIH5OIns75nFOe9ptjfpLVuOke+isY+NV4K4OBNSZmEr+kTwiygXEIO5JjA6he/3rqQ70SjBR1UhelnQLG/aWvddI68op4cBNtTCssmfEVNcWSWKZ+Q5xr9ccXvBIA5tLd/BLkFHfXpwHVnBPnfi+Rbj4Oqt/u3YXtmis3C5z5JbN68YelJIveEMCaU8gl5IaI/YRa06umgF3J3RpD/PU7vlKJkxQBy980JAC9rHCOf67VV+FY+SdDSk8Tq1fL9i5DiyDhoEXIi+isimiaid3WcLzW5rxYQlsGiCiYegdydERR8jJotU0TL3aT7Jr/rtbE8HVA/SQCc4lh0dFnkrwJ4StO51FgTQQrLYIkbTAzaX2EuNWy3tethUvK7Xltzu1VPDDdvFuAhlwlES2WnEOKnRDSo41xKwpYJLxQ1rKzyBBrT9+JWKaatamy27oH5Xa/75zgy4ojjwIAj4oX7M/URVCXKVZrFhoRqWZG4J3KE/A0hxBcVnw+jvoLBmjVrnjh9+nSs86/fvh2rrl9v2P5pfz8unToVe76m6esbx+DgUXR23sD8fB+uXt2N6emt2sfZsGEziBp/h0IQLl48t/R+/frtWLVK8vP7tB+XLqX7+c3OzqK7uzvVOXRj6nqLeK26GB/vw6FD6zA/v1yL0dGxgP37r2Dr1ukcZ2YO236fmzZtuiyEeLLhA1kDliQvONGkd6Psm6hplqnuQAbJpiFPRch/zBXffieFEG2+fdqEjsZTxWw8lLQxVzDFvFZ9LDe6Wixsoyud2Pb7hKJpVv5ZK1H93rZGkIwTpwrTX00ZVF1pO0Vb+MMO3DyBc+cu5trPnIlHvkIep3LC1giScaIK1giAO75td1DeYCdQrIU/GMYcutIP/xrAPwJYR0S/JqI/jHRgnMqJQuSHF5UogqWjqdcgipGDzjCMFy1CLoR4XgjxsHCcro8IIX4Y6cC4lRO554fXSZwGmUYM0wppmo6Ge+C0GihSDjrDMC75ulaK7PdWiXXiRhppClR0FLck7Wg4BuAHaGxydRvADgCD6OsbjzEPhmF0k6+Qm/Z7J7Wcg8Q6cSONNAUqQSvqRLXQkwb/RhDcqdBdCYmtc4bJi3yFXJffWybYaVrQBYl14kYaaXzUQSvqmHZ1hM8v/kpIDMPoJP/0w7R+b5Vg792bvAVdkFiHuYM8N5X127d7bhxpfNRR9gmz7qO4Z2R++KhuLg6aMkxe5C/kaVFZzzOKhRqitKALEusgd5DvprLq+nXPU0CaVXdkx8oIurYoPc1lQr8t4tiqni9e8mhkxTcOpvzYL+Rxe4NGCaQGiXWQOyjQf56mQMV/rGQ5O+fiAs4R5tpRCf0Z39hp/mTSxAmSCHIROz4yjH7sF3KVMPf2Jg+khvnu/e4gwHGnyDoOAZ6bTZoCFe+xxxHfug9z7QQJvXdsVeDzZsDY3nPF2e6KNyFZ+mPQjYMtdaY82C/kKuv5yJF0gdSovnuvO0VFrHTKKAKTxLoPc+1E9eEPAKNR9ouzj2y715oG5OmPYZa86gbh3gjYUmdKgqwBi+lXoqZZQSx3+hGZd/qpVOTNvNxXV1eM+SRp9HRSOA2yqP416b5Rxz7pXJdnv3v3OkLGjTuGEOpmYN5XWMM01TlaFdsroVdgW5OlpPB1FhMUtmmWDvKs+Azw0X/a3+88BQAR89nj+pDj+oCDXDtRrXz3/fJ+V67s92wPeqKI8yQRJfbhWvKqMVVPIQuK83HmDWMn5RDyPFG5TSqV5T7pkfPZ4/qQM1wFZ/Q5x0VF9Y6JNAWQAEa/4emzHuXGEjVOEOaqcd1CQWOqbhyVhGOiXsXKbhmmWLCQpyWsOlWVybJrl0TMVULSArn1p3MVnCBBHANG33A2u65qAUC0A6OPec6h88Yis6bdtrteSz5sTNmNI3kq6KOPHgsZj2Gyh4U8LWEZLirXy8KCxDJX5YurKjjTFBn5CRLEEQB3Jcf42+DqvLHIrOkTcH4Ok1i25JOMmTwVtKNDtVJOWZfMY2yAhVwHQT76oIyVhkrTKPniXusvTZGRnyBB9H12QHWczhsLEM0Nk3TMZKmg8/N9CcdjGHOUQ8jDmmMlbjurAZnrxUuDxe4VmEXVQZ59da2CEySIvs9GVcfpvLFEJdsxr17dnel4DBMF+4U8rDlW1OZZpsTedb20KqoxA3PMo1ibulbBCRLEGoA2yTHtWBawMSy7Z9xrzWJ5tWyXdHMCu7yEHFMwZDmJpl9a88hVedyVSrTPhXDyvLu6UuR/y1mRo5poDDMLCAePVxHqPPNezzx6lz5/772RjOeZH7blHSeFr7OYoLR55GFtZaO0nU3cYzwGiVr2Zr2AcFie+Q0sp67cWPo8OJMjKOea87EZRgf2C3lYW9koqxDF6TGexgWTqHCp+AsIqzM5gkrhuaEVw+jCfiEPy+OOsgqRKhjp355msYoSo87kaEVwSiPnYzOMDuwXcpXLAnAs5p07gc5OpxuiyqUxNyc/t397WhdMntkzBpFncgDBpfA6c84ZprmxX8gBeVtZr+U8M+OI8okTcpfGoiLNz7898TJvyNGaN++HXs7k6I14hCSlccVnDMPEoRxC7ieu5axKDfRvj+Jv1zUnLSRd3i0eTv+REQCKVZlW4E1p5HxshtFBOYU8ruU8PBxtexR/u645aSHp8m5xxHwM69YdwnLf8CBasZx1k3VGDsOUl3IKeVzL+eWXgRdeWLbAW1ud9y+/vHK/RCmEIWOvXq3Bb66yqpMu7xbnKWEEra3zEfddwHJKIlC8jJzwVMkNGzZLPmOYnJEll5t+JSoIirN4RJoCH42LVIQWBLW3C9HWlrIQKahoqCLkv4JK/VhSfB62YIMX1TmCXkUsFgr6OWZdmJU/thXKJMW264SiIKh4Qi4T0iTCHFWQvfv19jriqqnCs+GPxD+n3t6VY8mqTkOpCLVYhwlQ0LFpxw97xRkjLnFWTXKpCPU8gz4rJ7YJXFJsu047hFwl2FoEL+J4GscJ/SMhko9HOixi9xw6lncL4mR9qTf/Oe5XzMs/P90kvaagn6OOJxe7sE3gkmLbdaqEvFg+clVmx4wiGyJtoFA2nolxVKTJglneOWS7juXdgqjWl3rznyPs52oqzTCp3z9G98fQYxgmW7QIORE9RURXiOhDIvqTxCeKK5ixBC/FeGnHUZEmC2b5JJCn8W1DtLTC9AFHJ4/cf46gn5nJNMOkhUZh3R85VZIpLqmFnIhaAXwfwNMAHgfwPBE9nuhkKsHs7dUgeDHG0z2OijRZMMsnQaNVvQvAceTbx0S12lEvzKYZplloQvV0svyZEJwqyRQPHRb5VwB8KIS4KoS4A+AUgN9LdCaVhXrkiAbBizheW5tz4wD0jRNEokZaDSfBSov4DMLdC6YrPmXCeBLeronRiTPXNNZzmBtqEhcvnpN8xjD5Qo7/PMUJiH4fwFNCiN319zsB/CchxIu+/YbhmIVYs2bNE6dPn5aer298HI8eO4aO6WnM9/Xh6u7dmN66VbqvDrzjLXR14b5PPmnYZ3LXLkx+4xuxzz07O4vu7m4Ns4zHhg2bQdT4exWCcPHiOfT1jWPdukMr8r+FAO7e/Qw+/PCP666S6Ji8TtlcFxY6cOXKfuU8+/rG8eijx9DRMY35+T5cvbo79jWpyOt3mjV8ncVk06ZNl4UQTzZ8IIuAxnkB+H0AxzzvdwL430HHxE4/zANAvj3G/PKLiFdEcLqc6vMkWSumr7Migq8lW2zLckgKX2cxgcGslX8D8DnP+0fq2+ITpbFUFh0E3XP6x7CmjW2YeyEo8HcbwF4UZ8EH7pLIMGHoEPKfAXiMiD5PRO0AtgP4+0RnCmsslYWQumMAjWPk0vgqCWFphWGBvxkUZ8EHTv1jmDBSC7kQ4h6AFwH8GMD7AE4LId5LdLKwxlJZCGnQGFOKxlCq7bniBu5O1N/vxLJ1XYMj8FHJc8EHTv1jmDC05JELIc4IIdYKIf6jECL5f1hYgUwWHQSDxoja7lY3id1Jqu6GqL+PgylXRlhGCndJZJgwilXZGVYgo6USMoSgMRYUK96otusglTspqMqxojhG9SdhwpURtY1u0bokMkyxKJaQhxXIaKmElOC1eGdngfZ2+RgVhfj19poLwKZyJwUFClUui/+m2G7ClcHrdjKMDool5EBwgYyWSkgffot3Zsb5KlvjU1VA9PHH5gKwqdxJQYFClcviZcV2E1ZwGTJSzC+lxzBh3Jf3BGJTreqttJRZvHfvAt3dwI0bjWO7x1y75rhbZmcbm3q5FvOrr6af38CAPJgayZ1Ug+Oq8F6f17p2y8/9qLbrZgDylYVsyUhxXUPuz9cbg2D3D5MdxbPIsyauxet/Yrh5M97xcQlyJ4UGQbMOFMa1Tm3PSGHXEFMMWMjTBlBNB2BV7iQgYhBUR6AwikAnWf/T9oyUMriGmDLQnELuD262ta38PE4A1VQA1ossbpAoCJrEnxtVoJNapzZnpKQpVmLfOqOP5hNyWXCTSB7cjEKUAKzfBbJnz/L7hx5yXnEzXmIHQZNYzEB0gW5G6zSpayjp74Jh5NgX7EyKW2IvCxzeuSMPbkYlKADr3jhc63lqCnjlleXPvYFS1z3injOI2EHQIEEOGiuqQNseuEyC+3MbgfPzGIAj4mFGQNLfBcPIaQ6L3GuFqzC1nFvU5eRcouaIx3bpJLWYo7oPTAcui+qKSOIaasanF8YkzSHkUcTU1HJuSW4QUY6JnVOf1J8bVaBNBi7L5orgRmCMXppDyMOE0eRybkluEFGPibW6UFKLOY5Amwpcli3Nz/a0S6ZoNIeQBwmj6eXcZC6QIIzdVNJYzHlnlpTNFWF72iVTNJpDyFX+5JMnU6yT6UNVnFOtArt2Oe4PL+773t7kGTNx5uFMBnam+pXRFWHr74IpIs0h5CZ6tHhRdCjsGx93Pj9zxtnuRQhnHjduOK9Uiy8Hz0NfE6+8Ao7simCYIJpDyAFNq9UrUBTnPHrsmPN9Fn3UA+ahZ+GNPAOO7IpgmCCaR8hNohDkjulp55ss+qgHzEPPDSPvgCO7IhhGhb1CnsUizFFRCPJ8X5/zTRZl/AHz0HPDKFvAkWHKg51CXrTV7BVCfXX3bud70z76kHnouWGUMeDIMOXATiEv2mr2CqGe3rp15T6mfPQh89AzFgcc86eo1a1M3tgp5FkFD+PgFepaDRgZwYbNm7N3+xi7YXDAMV/KVt3K6MROIc8qeJgEj9uHiuD20QoHHPMj72AzU2TsFPKsgocqggKtUd0+RQrWMhbAwWZGjZ1CnlXwUEZYoDWK26dowVrGAjjYzKixU8iBbIKHMsIs7ihun6IFaxkL4GAzo8ZeIc+LMIs7ituniMFapuBwsJlRw0IelzCL2+P2ESq3T5GDtUyB4WAzI4eFPC5RLO662+fiuXNyt0/ewVqGYUpFcwi5zgwRHYHWsHNwRkvGcKENYzepFl8moj8AMArgtwB8RQjxlo5JaUW2+HHUBY5VBC22nPYcJubLBOAW2rjBZ7fQBmDXBWMLaS3ydwE8B+CnGuZiBtsyRGybr/VwoQ1jP6ksciHE+wBA/tVvioRtGSK2zdd6uNCGsZ9UQh4HIhrG8jPrPBG9m8W4XwKG2oB2//a7Qtx5h2jC8PAPAbgR54Cc55uU2NdZFL70JQy1tUl+3nfFnXfekf68rb3WmPB1FpOKbGOokBPROID/IPloRAjxo6ijCyGOwkl8BRG9JYR4MuqxtsLXWT6a5Vr5Ou0iVMiFEFvD9mEYhmHyoznSDxmGYUpMKiEnoq8T0a8B/GcAbxLRjyMeejTNuBbB11k+muVa+TotgoQQec+BYRiGSQG7VhiGYSyHhZxhGMZychNyIvoDInqPiBaJyPr0Hz9E9BQRXSGiD4noT/KejwmI6K+IaDqrmoC8IKLPEdF5IvrX+t/s3rznZAoiWkVE/0REP69f68G852QSImolon8hojfynksa8rTIi1/enxAiagXwfQBPA3gcwPNE9Hi+szLCqwCeynsSGXAPwEtCiMcBrAfwRyX9fQLAPIDNQojfBvBlAE8R0fp8p2SUvQDez3sSaclNyIUQ7wshruQ1vmG+AuBDIcRVIcQdAKcA/F7Oc9KOEOKnAG7mPQ/TCCH+XQjxz/XvP4bzj//ZfGdlBuEwW3/bVn+VMiOCiB4B8DUAx/KeS1rYR26GzwL4f573v0ZJ//GbDSIaBPA7AP5vzlMxRt3d8DaAaQA/EUKU9VoPA/gOnJU6rMaokBPROBG9K3mVzjplyg8RdQP4GwD7hBAf5T0fUwghFoQQXwbwCICvENEXc56SdojoGQDTQojLec9FB0abZjVxef+/Afic5/0j9W2MpRBRGxwRHxNCvJb3fLJACPEbIjoPJw5StoD2VwH8VyLaBmAVgM8Q0UkhxI6c55UIdq2Y4WcAHiOizxNRO4DtAP4+5zkxCSGnT/MPAbwvhPiLvOdjEiJaQ0QP1r/vBPC7AH6R66QMIIT4rhDiESHEIJz/z3O2ijiQb/ph0vL+wiOEuAfgRQA/hhMYOy2EeC/fWemHiP4awD8CWEdEvyaiP8x7Tob4KoCdADYT0dv117a8J2WIhwGcJ6J34BgkPxFCWJ2a1wxwiT7DMIzlsGuFYRjGcljIGYZhLIeFnGEYxnJYyBmGYSyHhZxhGMZyWMgZhmEsh4WcYRjGcv4/dYzmQVwfGDMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "visualisation(clusters, index, Pc, K=K)\n",
    "# print(Pc)\n",
    "# print(mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pc, index, img_seg = kmeans_image(path_image=path_image, K=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(103230,)\n",
      "[[252.69560225 251.67998184 249.86355505]\n",
      " [ 69.54698482  50.76095445  17.63930586]\n",
      " [228.4210873  174.61140875  69.5397857 ]\n",
      " [154.87751731  83.70097198  45.64258942]\n",
      " [         nan          nan          nan]]\n",
      "[[[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]\n",
      "\n",
      " [[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]\n",
      "\n",
      " [[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]\n",
      "\n",
      " ...\n",
      "\n",
      " [[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]\n",
      "\n",
      " [[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]\n",
      "\n",
      " [[252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  ...\n",
      "  [252 251 249]\n",
      "  [252 251 249]\n",
      "  [252 251 249]]]\n"
     ]
    }
   ],
   "source": [
    "print(index.shape)\n",
    "print(Pc)\n",
    "print(img_seg)"
   ]
  }
 ],
 "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
}