{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第9章 图像分割"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "平面设计师在对一些图像进行处理的时候，往往需要提取出图像中的某一个物体，如图 9-1 所示。这种技术被形象地称为“抠图”。“抠图”是图像分割的一个实际应用。在计算机视觉领域，图像分割是最古老也是被研究得最多的一个问题。根据是否基于人为给定的语义监督，图像分割可分为无监督图像分割和有监督图像分割。本章将聚焦于无监督图像分割，而有监督图像分割，包括语义分割和实例分割等，应归类于图像识别问题，将在后续章节介绍。无监督图像分割，其目的是把图像划分成若干互不相交的区域，使得具有相似特征的像素被划分到同一区域，不同特征的像素被划分到不同的区域（如图 9-2）。这一目标完全是通过挖掘图像内在结构信息实现，而不基于人为标注的语义类别，所以无监督图像分割是一个普适性的问题，其结果可以为解决高层视觉识别任务，如语义分割、目标检测、实例分割等，提供基础或者初始条件。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic2.zhimg.com/80/v2-451734557ce300fb4b053d7d59126335_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-1 “抠图”。</div>\n",
    "</center>\n",
    "\n",
    "在这一章中，我们将介绍两种具有代表性的无监督图像分割算法——基于K均值的图像分割算法和基于归一化图割的图像分割算法。\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-cc19aa6fad1171f69bbbf5a5896ec348_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-2 图像分割示例。</div>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.2 图像分割算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "既然本章聚焦于无监督图像分割，我们首先需要知道什么是无监督学习（unsupervised learning）。所谓无监督，指的是每一个用于训练的样本只有其本身的信息，而没有任何的标签信息，即没有“监督”信息。聚类（clustering）就是一种常用的无监督学习的方法，是指根据在数据中发掘出的样本与样本之间的关系，将数据进行分簇。其目的是使簇内的数据相互之间是相似的(相关的)，而不同簇中的样本是不同的(不相关的)。簇内相似性越大，簇间差距越大，说明聚类效果越好。作为一种无监督学习的方法，聚类是许多领域中常用的机器学习技术。图像分割本质是就是对像素进行聚类，所以图像分割一般都是基于聚类算法。接下来，我们将介绍两种代表性的基于聚类的图像分割算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.2.1 基于K均值聚类的图像分割算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "K均值（K-means）[1] [2] 是一个非常著名的聚类算法，它的基本思想是：\n",
    "\n",
    "1. 随机初始化$K$个聚类中心$c_1,...,c_K$；\n",
    "2. 给定聚类中心，确定每个聚类中的点：对于每个点$p$，我们计算与之距离最接近的$c_i$，并将$p$放入集群$i$中；\n",
    "3. 当所有点分配完毕时，将$c_i$设置为集群$i$中点的平均值；\n",
    "4. 如果有一个$c_i$发生了变化，请重复步骤2；反之，结束。\n",
    "\n",
    "同样，我们也可以将K均值算法应用于图像分割中，对图像中像素点进行聚类。那么，对于一张图像而言，我们该怎么定义像素点的信息呢？我们知道每个像素点包含有R、G、B三个通道的数值，如图 9-3 所示，因此我们可以以RGB值为特征对像素进行聚类。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-527eae0fc125d1987c9bf97e8e4a4f10_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-3 以RGB值为特征对像素进行聚类。</div>\n",
    "</center>\n",
    "\n",
    "我们来动手尝试编写kmeans算法，并利用RGB值对图像进行分割。我们先导入一张图像。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from matplotlib.image import imread\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from PIL import Image\n",
    "\n",
    "# 输入一张图像\n",
    "image = imread('t.jpeg')[:,:,:3]\n",
    "# 将RGB值统一到0-1内\n",
    "if np.max(image)>1:\n",
    "    image = image / 255\n",
    "\n",
    "X = image.reshape(-1, image.shape[2])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0 1 2 3]\n"
     ]
    }
   ],
   "source": [
    "# 利用kmeans算法进行聚类\n",
    "segmented_imgs = []\n",
    "\n",
    "# 设定聚类中心个数\n",
    "n_cluster= 4\n",
    "kmeans = KMeans(n_clusters=n_cluster, random_state=42).fit(X)\n",
    "print(np.unique(kmeans.labels_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-10-12T07:10:20.112994Z",
     "iopub.status.busy": "2022-10-12T07:10:20.111763Z",
     "iopub.status.idle": "2022-10-12T07:10:21.589547Z",
     "shell.execute_reply": "2022-10-12T07:10:21.588090Z",
     "shell.execute_reply.started": "2022-10-12T07:10:20.112945Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 159024 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43168 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 为每一个类别赋予一个对应的颜色，用于展示\n",
    "def decode_segmap(label_mask, plot=False):\n",
    "    label_colours = np.asarray([[79, 103, 67], [143, 146, 126], [129, 94, 64], [52, 53, 55],\n",
    "                       [96, 84, 70], [164, 149, 129]])\n",
    "    r = label_mask.copy()\n",
    "    g = label_mask.copy()\n",
    "    b = label_mask.copy()\n",
    "    # 为每个类别赋予对应的R、G、B值\n",
    "    for ll in range(0, 6):\n",
    "        r[label_mask == ll] = label_colours[ll, 0]\n",
    "        g[label_mask == ll] = label_colours[ll, 1]\n",
    "        b[label_mask == ll] = label_colours[ll, 2]\n",
    "\n",
    "    rgb = np.zeros((label_mask.shape[0], label_mask.shape[1], 3))\n",
    "    rgb[:, :, 0] = r\n",
    "    rgb[:, :, 1] = g\n",
    "    rgb[:, :, 2] = b\n",
    "\n",
    "    return rgb\n",
    "\n",
    "\n",
    "# 获得预测的标签\n",
    "segmented_img = kmeans.cluster_centers_[kmeans.labels_]\n",
    "segmented_imgs = decode_segmap(kmeans.labels_.reshape(image.shape[0],image.shape[1]))\n",
    "\n",
    "# 展示结果\n",
    "plt.imshow(image[:,:,:3])\n",
    "plt.title('Original image')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "\n",
    "plt.imshow(segmented_imgs.astype(np.uint8))\n",
    "plt.title('{} center'.format(n_cluster))\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们观察实验结果不难发现，仅仅利用RGB特征可能会在分割时出现一些离散的点。实际上，由于图像的像素点本身就蕴含着坐标信息，我们除了可以选择RGB值作为像素的特征，也可以将坐标信息纳入其中。而当我们在特征中纳入坐标信息后，我们便可以将图像中这些离散的点分割到周围的物体中。这种做法既考虑了颜色的相似性，也考虑了空间的离散程度。\n",
    "\n",
    "我们再次动手尝试图像分割，并比较设置不同个数的聚类中心对图像分割的影响。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-10-12T07:10:02.205430Z",
     "iopub.status.busy": "2022-10-12T07:10:02.204848Z",
     "iopub.status.idle": "2022-10-12T07:10:08.992311Z",
     "shell.execute_reply": "2022-10-12T07:10:08.991733Z",
     "shell.execute_reply.started": "2022-10-12T07:10:02.205401Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 242992 characters]\n",
      "text/plain": [
       "<Figure size 1200x800 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "image = imread('t.jpeg')[:,:,:3]\n",
    "# 将RGB值统一到0-255内\n",
    "if np.max(image)>1:\n",
    "    image = image / 255\n",
    "    \n",
    "\n",
    "sp = image.shape\n",
    "# 增加xy坐标的信息\n",
    "# 设定一个权重，对坐标信息加权\n",
    "weight = 2\n",
    "y = weight * np.array([[i for i in range(sp[1])] for j in range(sp[0])]) / sp[0] / sp[1]\n",
    "x = weight * np.array([[j for i in range(sp[1])] for j in range(sp[0])])/ sp[0] / sp[1]\n",
    "image = np.append(image, x.reshape(sp[0], sp[1], 1), axis=2)\n",
    "image = np.append(image, y.reshape(sp[0], sp[1], 1), axis=2)\n",
    "\n",
    "X = image.reshape(-1, image.shape[2])\n",
    "segmented_imgs = []\n",
    "\n",
    "# 将 K 分别设置为6, 5, 4, 3, 2\n",
    "n_colors = (6, 5, 4, 3, 2)\n",
    "for n_cluster in n_colors:\n",
    "    kmeans = KMeans(n_clusters=n_cluster, random_state=42).fit(X)\n",
    "    segmented_img = kmeans.cluster_centers_[kmeans.labels_]\n",
    "    segmented_imgs.append(decode_segmap(kmeans.labels_.reshape(image.shape[0],image.shape[1])).astype(np.uint8))\n",
    "\n",
    "# 展示结果\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.subplot(231)\n",
    "plt.imshow(image[:,:,:3])\n",
    "plt.title('Original image')\n",
    "plt.axis('off')\n",
    "\n",
    "for idx,n_clusters in enumerate(n_colors):\n",
    "    plt.subplot(232+idx)\n",
    "    plt.imshow(segmented_imgs[idx])\n",
    "    plt.title('{} center'.format(n_clusters))\n",
    "    plt.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们与之前没有加入坐标信息的特征分割的结果进行对比，可以很容易地发现加入了坐标信息之后，图像的分割变得更加平滑。除此之外，从实验结果中我们可以很明显的看出，随着$K$的增加，图像被分割的也越来越精细。实际上，$K$值的选择非常的重要，选择一个合适的$K$值可以帮助我们有效的进行图像分割。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.2.2 基于图割的算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于图割的方法利用图论领域的理论和方法，通常做法是将待分割的图像映射为带权无向图$\\mathbf{G}=(\\mathbf{V}，\\mathbf{E})$，其中，$\\mathbf{V}=\\{v_1, v_2,...,v_n\\}$是顶点的集合，$\\mathbf{E}$为边的集合。图中的每个节点$n \\in \\mathbf{V}$对应于图像中的每个像素，每条边$e \\in \\mathbf{E}$连接着一对相邻的像素，边的权值$w(v_i,v_j)$表示相邻像素之间在灰度、颜色或纹理方面的非负相似度，其中，$(v_i,v_j)\\in \\mathbf{E}$。因此，图像的分割问题便转化为图的切割（图割）问题，其中每一个子图对应着原图中被分割的区域，如图 9-5 所示。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic1.zhimg.com/80/v2-e05a7e6e1c42de7be9a33fa844c0ba04_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-5 图像分割与图的剪切。</div>\n",
    "</center>\n",
    "\n",
    "那么我们该如何对图进行切割呢？同kmeans算法的核心思想相同，我们希望划分后的子图内的顶点之间的相似度最大，而子图与子图之间的顶点对的相似度最小。我们以将一张图分割成两个子图为例，如图 9-6，把$\\mathbf{G}=(\\mathbf{V},\\mathbf{E})$的顶点分成两个子集$\\mathbf{A},\\mathbf{B}$；令$\\mathbf{A} \\cup \\mathbf{B} = \\mathbf{V}$，$\\mathbf{A} \\cap \\mathbf{B} = \\emptyset$，我们定义${\\rm cut}(\\mathbf{A},\\mathbf{B}) = \\sum_{\\mu \\in \\mathbf{A},v \\in \\mathbf{B}}w(\\mu, v)$ , 我们可以通过最小化${\\rm cut}(\\mathbf{A},\\mathbf{B})$实现上述思想。我们把这种算法称为最小割（MIN-CUT）[3]。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic4.zhimg.com/80/v2-4cf117da7a28862864eec6dcd12aab07_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-6 MIN-CUT分割图中节点。</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "然而当我们利用MIN-CUT分割图的时候，实际的分割情况往往如图 9-7 所示。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic4.zhimg.com/80/v2-9bbf1956fbd3e2e2fa8df8c2dc6ebbe7_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-7 MIN-CUT可能存在的问题。</div>\n",
    "</center>\n",
    "\n",
    "这种情况很容易解释，因为MIN-CUT是通过最小化${\\rm cut}(\\mathbf{A},\\mathbf{B})$实现分割，而当其中一个子图只存在一个节点的时候，${\\rm cut}(\\mathbf{A},\\mathbf{B})$很有可能是最小的。为了解决这种极端的分割的情况，Ncut（Normalized Cut）[4]算法应运而生。\n",
    "\n",
    "Ncut算法通过最小化下式对顶点进行分类：\n",
    "\n",
    "$${\\rm Ncut}(\\mathbf{A},\\mathbf{B}) = \\frac{{\\rm cut}(\\mathbf{A},\\mathbf{B})}{{\\rm assoc}(\\mathbf{A},\\mathbf{V})} + \\frac{{\\rm cut}(\\mathbf{A},\\mathbf{B})}{{\\rm assoc}(\\mathbf{B},\\mathbf{V})}$$\n",
    "\n",
    "其中，${\\rm assoc}(\\mathbf{A},\\mathbf{V})$表示$\\mathbf{A}$中所有顶点到图中所有点权重的和，即：\n",
    "\n",
    "$$ {\\rm assoc}(\\mathbf{A},\\mathbf{V}) = \\sum_{\\mu \\in \\mathbf{A},v \\in \\mathbf{V}}w(\\mu, v)$$\n",
    "\n",
    "我们不难发现，如果集合$\\mathbf{A}$只包含有一个点，那么${\\rm assoc}(\\mathbf{A},\\mathbf{V})$与${\\rm cut}(\\mathbf{A},\\mathbf{B})$是相等的，此时${\\rm Ncut}(\\mathbf{A},\\mathbf{B})$较大；而当集合$\\mathbf{A}$只包含有多个点时，那么${\\rm assoc}(\\mathbf{A},\\mathbf{V})$的值就会比${\\rm cut}(\\mathbf{A},\\mathbf{B})$要大，此时${\\rm Ncut}(\\mathbf{A},\\mathbf{B})$就会比较小。这也是为什么Ncut可以避免分割时MIN-CUT出现的问题。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"https://pic3.zhimg.com/80/v2-37f70f2b64c6a0217a805b5fe7aadb1a_1440w.jpg\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-8 向量表示顶点的聚类。</div>\n",
    "</center>\n",
    "\n",
    "为了更好的表示${\\rm Ncut}(\\mathbf{A},\\mathbf{B})$，如图 9-8 所示，我们用一个二值向量$\\mathbf{x}\\in \\{-1,+1\\}^n$表示顶点的分布情况：若顶点处于黑色位置，则对应的值为1，反之为-1。如此，我们可以将${\\rm Ncut}(\\mathbf{A},\\mathbf{B})$改写为：\n",
    "\n",
    "$${\\rm Ncut}(\\mathbf{A},\\mathbf{B}) = \\frac{\\sum_{(x_{i}>0, x_{j}<0)}-w_{i j} x_{i} x_{j}}{\\sum_{x_{i}>0} d_{i}}+\\frac{\\sum_{(x_{i}<0, x_{j}>0)}-w_{i j} x_{i} x_{j}}{\\sum_{x_{i}<0} d_{i}}$$\n",
    "\n",
    "其中，$d_i = \\sum_j w(i,j)$。我们令$\\mathbf{D}$为对角矩阵，其中$\\mathbf{D}_{ii} = d_i$，令\n",
    "$\\mathbf{W}$为权重矩阵，即$\\mathbf{W}_{ij} = w(i,j)$。令\n",
    "\n",
    "$$\\mathbf{y} = \\frac{1}{2}[(\\mathbf{1}+\\mathbf{x}) - b(\\mathbf{1}-\\mathbf{x})] \\in \\{1, -b\\}^n$$\n",
    "\n",
    "其中$b=\\frac{\\sum_{x_i>0}d_i}{\\sum_{x_i<0}d_i}$，我们有\n",
    "\n",
    "$$\\min_{\\mathbf{x}} {\\rm Ncut}(\\mathbf{x}) = \\min_{\\mathbf{y}} \\frac{\\mathbf{y}^\\mathsf{T}(\\mathbf{D}-\\mathbf{W})\\mathbf{y}}{\\mathbf{y}^\\mathsf{T}\\mathbf{D}\\mathbf{y}}$$\n",
    "\n",
    "上述优化问题的最优解可以通过下式求得：\n",
    "\n",
    "$$(\\mathbf{D}-\\mathbf{W})\\mathbf{y} = \\lambda \\mathbf{D}\\mathbf{y}$$\n",
    "\n",
    "即：\n",
    "\n",
    "$$\\mathbf{D}^{-\\frac{1}{2}}(\\mathbf{D}-\\mathbf{W})\\mathbf{D}^{{-\\frac{1}{2}}}\\mathbf{D}^{{\\frac{1}{2}}}\\mathbf{y} = \\lambda \\mathbf{D}^{{\\frac{1}{2}}}\\mathbf{y}$$\n",
    "\n",
    "令$\\mathbf{z}=\\mathbf{D}^{{\\frac{1}{2}}}\\mathbf{y}$，则：\n",
    "\n",
    "$$\\mathbf{D}^{-\\frac{1}{2}}(\\mathbf{D}-\\mathbf{W})\\mathbf{D}^{{-\\frac{1}{2}}}\\mathbf{z} = \\lambda\\mathbf{z}$$\n",
    "\n",
    "很明显，这里$\\mathbf{z}$和$\\lambda$分别是$\\mathbf{D}^{-\\frac{1}{2}}(\\mathbf{D}-\\mathbf{W})\\mathbf{D}^{{-\\frac{1}{2}}}$的特征向量和特征值。因此，我们只需求解该特征值即可。不过需要注意的是，我们要的分割点是第二小的特征值所对应的特征向量，因为以最小的特征值对应的特征向量作为划分依据则会将图中所有的顶点聚为一类，没有任何的意义。\n",
    "\n",
    "我们把Ncut算法的计算过程总结为以下几步：\n",
    "\n",
    "1. 输入图像，并以每个像素点作为一个节点构建图；\n",
    "\n",
    "2. 计算权重矩阵$\\mathbf{W}$和对角矩阵$\\mathbf{D}$，这里我们将$\\mathbf{W}$定义为：\n",
    "$$\\begin{aligned} w_{i j}=& \\exp({\\frac{-\\left\\|\\mathbf{I}_{i}-\\mathbf{I}_{j}\\right\\|_{2}^{2}}{\\sigma_{I}^{2}} }) * \\begin{cases}\\exp({\\frac{-\\left\\|{x}_i-{x}_j\\right\\|_{2}^{2}}{\\sigma_{X}^{2}}}) & \\text { if }\\|x_i-x_j\\|_{2}<r \\\\ 0 & \\text { otherwise. }\\end{cases} \\end{aligned}$$\n",
    "其中，$\\sigma_{\\mathbf{I}}$，$\\sigma_{\\mathbf{X}}$和$r$是三个超参数；\n",
    "\n",
    "3. 计算$\\mathbf{D}^{-\\frac{1}{2}}(\\mathbf{D}-\\mathbf{W})\\mathbf{D}^{{-\\frac{1}{2}}}$的特征值和特征向量，并确定第二小的特征值所对应的特征向量。\n",
    "\n",
    "接下来，我们将动手编程实现Ncut算法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-20T06:59:29.560525Z",
     "iopub.status.busy": "2022-11-20T06:59:29.559668Z",
     "iopub.status.idle": "2022-11-20T06:59:30.809302Z",
     "shell.execute_reply": "2022-11-20T06:59:30.808285Z",
     "shell.execute_reply.started": "2022-11-20T06:59:29.560394Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ab257610>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 98900 characters]\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sample_img = cv2.imread('t.jpeg', cv2.IMREAD_GRAYSCALE).astype(np.float32) / 255\n",
    "plt.imshow(sample_img, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们接着编写权重矩阵。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-20T06:59:30.812754Z",
     "iopub.status.busy": "2022-11-20T06:59:30.811390Z",
     "iopub.status.idle": "2022-11-20T07:24:25.257332Z",
     "shell.execute_reply": "2022-11-20T07:24:25.256304Z",
     "shell.execute_reply.started": "2022-11-20T06:59:30.812708Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.81873075, 0.44863849, ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.81873075, 1.        , 0.81747264, ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       [0.44863849, 0.81747264, 1.        , ..., 0.        , 0.        ,\n",
       "        0.        ],\n",
       "       ...,\n",
       "       [0.        , 0.        , 0.        , ..., 1.        , 0.81873075,\n",
       "        0.44863849],\n",
       "       [0.        , 0.        , 0.        , ..., 0.81873075, 1.        ,\n",
       "        0.81747264],\n",
       "       [0.        , 0.        , 0.        , ..., 0.44863849, 0.81747264,\n",
       "        1.        ]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def cal_dist_weighted_matrix(size, r, sig_X):\n",
    "    h, w = size\n",
    "    X_matrix = np.zeros((h*w, h*w))\n",
    "    for i in range(h*w):\n",
    "        for j in range(h*w):\n",
    "            i_row, i_col = i // w, i % w\n",
    "            j_row, j_col = j // w, j % w\n",
    "            dist = np.power(i_row - j_row, 2) + np.power(i_col - j_col, 2)\n",
    "            \n",
    "            if np.sqrt(dist) < r:    \n",
    "                X_matrix[i, j] = np.exp(-dist / sig_X)\n",
    "    return X_matrix\n",
    "    \n",
    "    \n",
    "def set_weighted_matrix(img, sig_I=0.01, sig_X=5, r=10):\n",
    "    vec_img = img.flatten()\n",
    "    F_matrix = np.power(vec_img[None, :] - vec_img[:, None], 2)\n",
    "    F_matrix /= sig_I\n",
    "    F_matrix = np.exp(-F_matrix)\n",
    "    \n",
    "    X_matrix = cal_dist_weighted_matrix(img.shape, r, sig_X)\n",
    "    return F_matrix * X_matrix\n",
    "\n",
    "\n",
    "weighted_matrix = set_weighted_matrix(sample_img)\n",
    "\n",
    "weighted_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们接着求解第二大的特征值对应的特征向量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-20T07:29:40.250671Z",
     "iopub.status.busy": "2022-11-20T07:29:40.250227Z",
     "iopub.status.idle": "2022-11-20T07:29:40.259628Z",
     "shell.execute_reply": "2022-11-20T07:29:40.258310Z",
     "shell.execute_reply.started": "2022-11-20T07:29:40.250634Z"
    }
   },
   "outputs": [],
   "source": [
    "def n_cuts(W, image):\n",
    "    n = W.shape[0]\n",
    "    s_D = sparse.csr_matrix((n, n))\n",
    "    d_i = W.sum(axis=0)\n",
    "    for i in range(n):\n",
    "        s_D[i, i] = d_i[i]\n",
    "    s_W = sparse.csr_matrix(W)\n",
    "    print(s_W.shape)\n",
    "    s_D_nhalf = np.sqrt(s_D).power(-1)\n",
    "    print(s_D.shape)\n",
    "\n",
    "    L = s_D_nhalf @ (s_D - s_W) @ s_D_nhalf\n",
    "    print(L.shape)\n",
    "\n",
    "    _, eigenvalues, eigenvectors = svds(L, which='SM')\n",
    "    print(eigenvectors.shape)\n",
    "\n",
    "    for i in range(1, 8):\n",
    "        print(eigenvectors[i].shape)\n",
    "        Partition = eigenvectors[i] > np.sum(eigenvectors[i])/len(eigenvectors[i])\n",
    "        print(Partition.shape)\n",
    "        skimage.io.imshow(Partition.reshape(image.shape))\n",
    "        plt.title('Ncut')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们展示相应的分割结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(17160, 17160)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weighted_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-11-20T07:29:43.523582Z",
     "iopub.status.busy": "2022-11-20T07:29:43.523143Z"
    }
   },
   "outputs": [],
   "source": [
    "from scipy import sparse\n",
    "from scipy.sparse.linalg import eigs, svds\n",
    "import skimage\n",
    "n_cuts(weighted_matrix, sample_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import skimage\n",
    "import skimage.io\n",
    "from scipy import sparse\n",
    "from scipy.sparse.linalg import eigs,svds\n",
    "import matplotlib.pyplot as plt\n",
    "from time import *\n",
    "\n",
    "def generate_W(img):\n",
    "    N = len(img)\n",
    "    d = np.zeros((N, N), dtype=np.uint8)\n",
    "    for i in range(N):\n",
    "        d[i, :] = np.abs(img - img[i])\n",
    "    W = np.exp(-np.power(d / np.max(d), 2))\n",
    "    return W\n",
    "\n",
    "def n_cuts(W):\n",
    "    D = np.sum(W, axis=1)\n",
    "    D = np.diag(D)\n",
    "\n",
    "    s_D = sparse.csr_matrix(D)\n",
    "    s_W = sparse.csr_matrix(W)\n",
    "    s_D_nhalf = np.sqrt(s_D).power(-1)\n",
    "    L = s_D_nhalf @ (s_D - s_W) @ s_D_nhalf\n",
    "    _, eigenvalues, eigenvectors = svds(L,which='SM')#svd much faster than eig(lanczos method)\n",
    "\n",
    "    plt.figure()\n",
    "    for i in range(5):\n",
    "        Partition = eigenvectors[i] > 0#choose 0 as split point, which can also be replaced by np.median(eigenvectors[i+1])\n",
    "        plt.subplot(1, 5, i+1)\n",
    "        skimage.io.imshow(Partition.reshape(image.shape))\n",
    "    plt.show()\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    image = skimage.img_as_ubyte(skimage.io.imread('t.jpeg',as_gray=True))#smaller than 100*100px better\n",
    "    W = generate_W(image.flatten())\n",
    "    n_cuts(W)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.3 小结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本章中我们学习了图像分割的相关方法。图像分割是计算机视觉中的一个重要操作，也是实例分割、语义分割的基础，在之后的相关章节中我们也会进行相应的介绍。至此，基础的图像处理我们已经介绍完毕，在接下来的章节中，我们会把目光聚焦在图像的语义理解上。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.4 习题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于Ncut算法，你还可以自己设计出更加有效的权重矩阵$\\mathbf{W}$吗？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] A. K. Jain and R. C. Dubes, Algorithms for Clustering Data, Englewood Cliffs, NJ, USA:Prentice-Hall, 1988.\n",
    "\n",
    "[2] J. MacQueen, \"Some methods for classification and analysis of multivariate observations\", Proc. 5th Berkeley Symp. Math. Statist. Probab., vol. 1, pp. 281-297, 1967.\n",
    "\n",
    "[3] Mechthild Stoer, Frank Wagner: A simple min-cut algorithm. J. ACM 44(4): 585-591 (1997)\n",
    "\n",
    "[4] Jianbo Shi, Jitendra Malik: Normalized Cuts and Image Segmentation. IEEE Trans. Pattern Anal. Mach. Intell. 22(8): 888-905 (2000)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
