{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第9章 图像分割"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.1 简介"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "所谓图像分割（Image Segmentation），是指根据灰度、彩色、空间纹理、几何形状等特征把图像划分成若干个互不相交的区域，使得这些特征在同一区域内表现出一致性或相似性，而在不同区域间表现出明显的不同。简单的说就是在一副图像中，把目标从背景中分离出来，如图 9-1 所示。对于灰度图像来说，区域内部的像素一般具有灰度相似性，而在区域的边界上一般具有灰度不连续性。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"1.png\" width=500>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-1 图像分割示例</div>\n",
    "</center>\n",
    "\n",
    "图像分割在计算机视觉中有着极为重要的作用。今天，就让我们一起来学习图像分割的相关方法吧。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.2 图像分割算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "图像分割的算法步骤主要分为两步：\n",
    "\n",
    "1. 自底向上：将具有相似特征的标记分组；\n",
    "2. 自顶向下：对可能属于同一对象的标记进行分组。\n",
    "\n",
    "我们将介绍两种不同的分割算法，分别基于聚类和基于图。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.2.1 基于聚类的算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "聚类是根据在数据中发现的描述对象及其关系的信息，将数据对象分簇。其目的是使簇内的对象相互之间是相似的(相关的)，而不同簇中的对象是不同的(不相关的)。簇内相似性越大，簇间差距越大，说明聚类效果越好。\n",
    "\n",
    "聚类是一种无监督学习的方法，是许多领域中常用的机器学习技术，而K-means是一个非常著名的聚类算法，它的基本思想是：\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-means算法的图像分割。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 203672 characters]\n",
      "text/plain": [
       "<Figure size 864x576 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from matplotlib.image import imread\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "image = imread('graf.png', 0)\n",
    "X = image.reshape(-1, 3)\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(segmented_img.reshape(image.shape))\n",
    "\n",
    "# 展示结果\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.subplot(231)\n",
    "plt.imshow(image.astype('uint8'))\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].astype('uint8'))\n",
    "    plt.title('{} colors'.format(n_clusters))\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从实验结果中我们可以很明显的看出，随着$K$的增加，图像被分割的也越来越精细。实际上，$K$值的选择非常的重要，选择一个合适的$K$值可以帮助我们有效的进行图像分割。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9.2.2 基于图的算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于图论的方法利用图论领域的理论和方法，将图像映射为带权无向图，因此我们便将图像分割问题转化为图的顶点划分问题，之后便可以利用最小剪切准则得到图像的最佳分割。此类方法把图像分割问题与图的最小割（MIN-CUT）问题相关联，通常做法是将待分割的图像映射为带权无向图$\\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-2 所示。\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"2.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-2 图像分割与图的剪切</div>\n",
    "</center>\n",
    "\n",
    "分割的原则就是使划分后的子图在内部保持相似度最大，而子图之间的相似度保持最小。我们以一个两类的分割为例，把$\\mathbf{G}=(\\mathbf{V}，\\mathbf{E})$分成两个子集$\\mathbf{A},\\mathbf{B}$；令$\\mathbf{A} \\cup \\mathbf{B} = \\mathbf{V}$，$cut(\\mathbf{A},\\mathbf{B}) = \\sum_{\\mu \\in \\mathbf{A},v \\in \\mathbf{B}}w(\\mu, v)$ , 其中$w(\\mu, v)$是权重, 最小割就是让$cut(\\mathbf{A},\\mathbf{B})$值最小的分割。Ncut（Normalized Cut）算法是一个代表性的基于图的图像分割算法。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"3.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-3 MIN-CUT分割图中节点</div>\n",
    "</center>\n",
    "\n",
    "\n",
    "要想理解Ncut与MIN-CUT之间的关系，我们先看图 9-3。我们假设当满足最小割时，红色的边就是分割的边界，但当我们实际运用最小割的时候，结果往往如图 9-4 所示。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"4.png\" width=400>\n",
    "    <br>\n",
    "    <div style=\"color:orange; \n",
    "    display: inline-block;\n",
    "    color: #999;\n",
    "    padding: 2px;\">图9-4 MIN-CUT可能存在的问题</div>\n",
    "</center>\n",
    "\n",
    "这种情况很容易理解，因为最小化切割就是让$cut\\mathbf{A},\\mathbf{B})$的值最小的情况，而边缘处其值确实是最小。为了解决这个情况，我们便需要应用Ncut算法。\n",
    "\n",
    "Ncut算法通过最小化下式对节点进行分类：\n",
    "\n",
    "$$Ncut(\\mathbf{A},\\mathbf{B}) = \\frac{cut(\\mathbf{A},\\mathbf{B})}{assoc(\\mathbf{A},\\mathbf{V})} + \\frac{cut(\\mathbf{A},\\mathbf{B})}{assoc(\\mathbf{B},\\mathbf{V})}$$\n",
    "\n",
    "其中，$assoc(\\mathbf{B},\\mathbf{V})$表示$\\mathbf{A}$中所有点到图中所有点的权重的和，即：\n",
    "\n",
    "$$ assoc(\\mathbf{A},\\mathbf{V}) = \\sum_{\\mu \\in \\mathbf{A},v \\in \\mathbf{V}}w(\\mu, v)$$\n",
    "\n",
    "我们不难发现，如果集合$\\mathbf{A}$只包含有一个点，那么$assoc(\\mathbf{A},\\mathbf{V})$与$cut(\\mathbf{A},\\mathbf{B})$是相等的；而当集合$\\mathbf{A}$只包含有多个点时，那么$assoc(\\mathbf{A},\\mathbf{V})$就会比$cut(\\mathbf{A},\\mathbf{B})$要大。这也是为什么Ncut可以避免分割时出现一个点的情况。\n",
    "\n",
    "\n",
    "<center>\n",
    "    <img style=\"border-radius: 0.3125em;\" \n",
    "    src=\"5.png\" 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",
    "那么我们该如何用向量表示聚类呢？如图 9-5 所示，我们用一个二值向量$x\\in \\{-1,+1\\}^n$表示聚类情况：若节点处于黑色部分，则对应的值为1，反之为0。因此，我们可以将$Ncut(\\mathbf{A},\\mathbf{B})$改写为：\n",
    "\n",
    "$$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}} 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算法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.linalg.decomp import eig\n",
    "from scipy import sparse\n",
    "from scipy.sparse.linalg import eigs\n",
    " \n",
    "# 定义Ncut类\n",
    "class Ncut(object):\n",
    "    def __init__(self, img):\n",
    "        # 获得输入图像的基本信息\n",
    "        self.no_rows, self.no_cols, self.channel = img.shape\n",
    "        # 获得节点的总数\n",
    "        self.N = self.no_rows * self.no_cols\n",
    "        self.V_nodes = self.V_node_maker(img)\n",
    "        self.X = self.X_maker()\n",
    "        self.F = self.F_maker(img)\n",
    "        # parameter for W clculate\n",
    "        self.r = 2\n",
    "        self.sigma_I = 4\n",
    "        self.sigma_X = 6\n",
    "        # Dense W,D\n",
    "        self.W = self.W_maker()\n",
    "        self.D = self.D_maker()\n",
    " \n",
    "    # V_nodes shape : [self.N, 1, 3]\n",
    "    def V_node_maker(self, img):\n",
    "        b, g, r = cv2.split(img)\n",
    "        b = b.flatten()\n",
    "        g = g.flatten()\n",
    "        r = r.flatten()\n",
    "        V_nodes = np.vstack((b,g))\n",
    "        V_nodes = np.vstack((V_nodes, r))\n",
    "        return V_nodes\n",
    " \n",
    "    def X_maker(self):\n",
    "        X_temp = np.arange(self.N)\n",
    "        X_temp = X_temp.reshape((self.no_rows,self.no_cols))\n",
    "        X_temp_rows = X_temp // self.no_rows\n",
    "        X_temp_cols = (X_temp // self.no_cols).T\n",
    "        X = np.zeros((self.N, 1, 2))\n",
    "        X[:,:,0] = X_temp_rows.reshape(self.N,1)\n",
    "        X[:,:,1] = X_temp_cols.reshape(self.N,1)\n",
    "        return X\n",
    " \n",
    "    def F_maker(self,img):\n",
    "        if self.channel < 2:\n",
    "            return self.gray_feature_maker(img)\n",
    "        else:\n",
    "            return self.color_img_feature_maker(img)\n",
    " \n",
    "    def gray_feature_maker(self,img):\n",
    "        print('need to ')\n",
    " \n",
    "    def color_img_feature_maker(self,img):\n",
    "        F = img.flatten().reshape((self.N,1,self.channel))\n",
    "        F = F.astype('uint8')\n",
    "        return F\n",
    " \n",
    "    def W_maker(self):\n",
    "        X = self.X.repeat(self.N,axis = 1)\n",
    "        X_T = self.X.reshape((1,self.N,2)).repeat(self.N,axis = 0)\n",
    "        diff_X = X - X_T\n",
    "        diff_X = diff_X[:,:,0]**2 + diff_X[:,:,1]**2\n",
    " \n",
    "        F = self.F.repeat(self.N,axis = 1)\n",
    "        F_T = self.F.reshape((1,self.N,3)).repeat(self.N,axis = 0)\n",
    "        diff_F = F - F_T\n",
    "        diff_F = diff_F[:,:,0]**2 + diff_F[:,:,1]**2 + diff_F[:,:,2]**2\n",
    " \n",
    "        W_map = diff_X < self.r**2 # valid map for W\n",
    " \n",
    "        W = np.exp(-((diff_F / (self.sigma_I**2)) + (diff_X / (self.sigma_X**2))))\n",
    " \n",
    "        return W * W_map \n",
    " \n",
    "    def D_maker(self):\n",
    "        # D is a diagonal matrix using di as diagonal, di is the sum of weight of node i with all other nodes\n",
    "        d_i = np.sum(self.W, axis=1)\n",
    "        D = np.diag(d_i)\n",
    "        return D\n",
    " \n",
    "    def EigenSolver(self):\n",
    "        L = self.D - self.W\n",
    "        R = self.D\n",
    "        lam,y = eig(L, R)\n",
    "        index = np.argsort(lam)\n",
    " \n",
    "        top2 = lam[index[:2]].real\n",
    "        smallest_2 = y[:,index[1]]\n",
    "        print('dense eigenvector')\n",
    "        return smallest_2.real\n",
    " \n",
    "    def EigenSolver_sparse(self):\n",
    "        s_D = sparse.csr_matrix(self.D)\n",
    "        s_W = sparse.csr_matrix(self.W)\n",
    "        s_D_nhalf = np.sqrt(s_D).power(-1)\n",
    "        L = s_D_nhalf @ (s_D - s_W) @ s_D_nhalf\n",
    "        lam,y = eigs(L)\n",
    "        index = np.argsort(lam)\n",
    " \n",
    "        top2 = lam[index[:2]].real\n",
    "        smallest_2 = y[:, index[1]]\n",
    "        print('sparse eigenvector')\n",
    "        return smallest_2.real\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread('graf.png', cv2.IMREAD_COLOR)\n",
    "cutter = Ncut(img)\n",
    "eigenvector = cutter.EigenSolver_sparse()\n",
    "b = eigenvector\n",
    "b = b.reshape((50,50)).astype('float64')\n",
    "b = (b / b.max())*255\n",
    "plt.imshow(b.astype('uint8'))\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.3 总结"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在本章中我们学习了图像分割的相关方法。图像分割是计算机视觉中的一个重要操作，也是实例分割、语义分割的基础，在之后的相关章节中我们也会进行相应的介绍。至此，基础的图像处理我们已经介绍完毕，接下来，我们会把目光聚焦在图像的语义理解。"
   ]
  }
 ],
 "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": 2
}
