图割算法之Graph Cuts浅见

前言

　　以后可能需要用图割算法做一些事情，所以就简单阅读了一下图割算法中最基础的一篇论文，Boykov Y Y, Jolly M P. Interactive graph cuts for optimal boundary & region segmentation of objects in ND images. Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on. IEEE, 2001, 1: 105-112. 。

图篇

　　图割算法（Graph Cuts），顾名思义，是基于图的一种图像分割方法，这里的图既是指数据结构中的图结构也是指数学中的图论，毕竟数据结构中的图结构也来自于数学中图论。既然是图，必然离不开图的构建，论文中图的构建如下，可以简称为 S/T 图：

和普通的无向图 \(G=(V,E)\) 一样，该图也是由顶点集合 \(V\) 和边集合 \(E\) 构成，不一样的是，该图有两种顶点和两种边：

第一种顶点是普通顶点。由图像中各像素点组成，每个顶点对应于图像中每个像素点。每两个邻域顶点，对应于图像中每相邻两个像素点（对2维图像来说，是8邻域；对3维图像来说，是26邻域），的连接就是一条边，这种边叫做 n-links (neighborhood links)。
第二种顶点是两个终端顶点，代表目标（object）的 source terminal（简称 S）和代表背景（background）的 sink terminal（简称 T），每个终端顶点都与所有普通顶点相连，这种相连构成的边叫做 t-links (terminal links)。

　　设图像 \(P\) 中每个像素点为 \(p\)，则 \(p \in P\)，目标终端顶点 \(S\)，背景终端顶点 \(T\) ，则构建的图 \(G\) 中顶点集合可表示为：\(V = P \cup \{S,T\}\) ；每个 \(p\) 都有两种 t-links 边，分别为 \(\{p,S\}\) 和 \(\{p,T\}\) ，每个 \(p\) 与其邻域内像素 \(q\) 构成的 n-links 边为 \(\{p,q\}\) ，设 \(N\) 代表 n-links 边集合，即邻域边集合，则 \(\{p,q\} \in N\) ，则构建的图 \(G\) 中边集合可表示为：\(E = N \bigcup\limits_{p \in P} \{\{p,S\},\{p,T\}\}\) 。

因为这构建的是一个无向带权图，所以还需要设定边的权值，论文中边的权值设定如下：

边	权值（代价，能量）	区域
\(\{p,q\}\)	\(B_{\{p,q\}}\)	\(\{p,q\} \in N\)
\(\{p,S\}\)	\(\lambda \cdot R_p("bkg")\)	\(p \in P, p \notin O \cup B\)
\(\{p,S\}\)	\(K\)	\(p \in O\)
\(\{p,S\}\)	0	\(p \in B\)
\(\{p,T\}\)	\(\lambda \cdot R_p("obj")\)	\(p \in P, p \notin O \cup B\)
\(\{p,T\}\)	0	\(p \in O\)
\(\{p,T\}\)	\(K\)	\(p \in B\)

其中集合 \(O\) 代表确定为目标的像素点的集合，即在交互时手动标注为目标的像素点；
集合 \(B\) 代表确定为背景的像素点的集合，即在交互时手动标注为背景的像素点；
\(B_{\{p,q\}}\) 可以表示 {p,q} 不连续的惩罚因子，其正比于一个指数函数，具体为：\(B_{\{p,q\}} \propto exp\left(-\frac{(I_p-I_q)^2}{2\sigma^2}\right) \cdot \frac{1}{dist(p,q)}\) ，其中 \(I_p\) 和 \(I_q\) 代表像素 \(p\) 和像素 \(q\) 的像素值，\(dist(p,q)\) 代表其两像素点之间的距离。从 \(B_{\{p,q\}}\) 对应的函数中可以看出当两像素点之间的差异越大时，其权值越小；两像素点之间距离越大时，权值越小（这里对于二维图像来说，一般只考虑周围 8 邻域内像素点，因此可以简单认为相邻两像素点之间距离为 1）。
\(R_p("bkg")\) 和 \(R_p("obj")\) 分别是指像素 p 分配给背景和前景目标的惩罚因子，该惩罚因子与像素 p 属于前景目标的概率 \(Pr(I|O)\) 和背景的概率 \(Pr(I|B)\) 有关，一般取概率的负对数值，具体如下：\(R_p(“obj”) = −lnPr(I_p|O)\) 和 \(R_p(“bkg”) = −lnPr(I_p|B)\) 。至于概率的计算方法可以用简单的直方图概率模型，因为 \(O\) 和 \(B\) 为手动标注的确定的前景目标和背景的像素点集合，因此可以分别统计其灰度直方图，再对直方图频数进行归一化得到分布概率直方图，将像素 p 与背景分布概率直方图和前景目标概率分布直方图进行对比，即可得到其属于背景的概率和属于前景目标的概率（当然更好的一种计算前景和背景概率的方法是用高斯混合模型）。
\(K\) 表示整个无向带权图中最大的权值，具体计算方法如下：\(K=1+\max\limits_{p \in P} \sum \limits_{q:\{p,q\} \in N} B_{\{p,q\}}\) ，即取图像中像素点邻域边权值之和的最大值再加 1 （对于二维图像而言，邻域边权值之和是指 8 邻域边权值之和，图像中每个像素点都有其邻域边权值之和，取所有像素点中的最大值），至于为什么要令 K 为整个无向带权图中权值最大值？请详见下文。
至于其中的参数 \(\lambda\) 和 \(\sigma\) 论文中没有明确指定，在初始计算时可以将其置为 1，随后再慢慢调整。

至此，整个 Graph Cuts 算法中关于图的构建已经完全确定，接下来就是割（cut）的算法了。

割篇

　　割（cut）是构建的无向带权图中边集合 \(E\) 的一个子集 \(C\)，即割 \(C \subset E\) ，该 cut 的代价（cost）可表示为：\(|C|=\sum \limits_{e \in C} W_e\)，即该割边集合的所有边权值之和。

如果一个割，其包含的所有边的权值之和最小，那么这个割就称为最小割，也就是图割的结果。在图分割里，是要求两个终端顶点被分离开的，即最小割把图的顶点划分为两个不相交的子集 \(S\) 和 \(T\) ，其中 \(s∈S\)，\(t∈T\)，\(T=V/S\)。事实上，这两个子集对应于图像的前景像素集和背景像素集，也就相当于完成了图像分割。

　　割相关的算法有很多，eg：Max-flow/Min-cut、GrabCut、One cut、 Normalized cut、Ratio cut 等等，但是其最关键的地方在于能量方程（或称代价函数，损失函数）的优化，能量优化的目的在于最小化能量函数，而最小化能量函数时取得割就是最小割，即图割的结果。论文中能量方程的定义如下： \[ E(A)= \lambda \cdot R(A)+B(A) \] 其中： \[ R(A)= \sum_{p \in P} R_p(A_p) \\ B(A)= \sum_{\{p,q\} \in N} B_{\{p,q\}} \cdot \delta(A_p,A_q) \\ 其中　\delta(A_p,A_q) = \begin{cases} 1 & \text{if } A_p \neq A_q \\ 0 & \text{otherwise}. \end{cases} \]

其中令 \(A=(A_1,\cdots,A_p,\cdots,A_{|P|})\) 为二值向量，每个 \(A_p\) 可赋值为 “obj”（可用 1 表示前景）或 “bkg”（可用 0 表示背景）；
\(R(A)\) 表示区域项，即上文中的两种 t-links 边相连的图像像素点权值，代表像素点分配给 “obj” 或 “bkg” 的惩罚项，对应于 \(R_p("obj")\) 和 \(R_p("bkg")\) ，求和后即可得到 \(R(A)\) 。当像素点 p 属于前景目标的概率 \(Pr(I|O)\) 大于背景的概率 \(Pr(I|B)\) 时，由上文 \(R_p(A_p) = −lnPr(I_p|A_p)\) 可知，此时 \(R_p("obj")\) 小于 \(R_p("bkg")\) ，即当像素 p 更有可能属于目标时，将 p 归类为目标就会使能量 \(R(A)\) 小，也因此要取概率的负对数值，由此，如果所有像素点都能正确划分为前景和背景，则总能量会达到最小。当像素点 p 已经人为手动标注为确定的前景目标时，这时其与顶点 \(S\) 相连的边的权值为 K，为整个无向带权图中最大的权值，即能量也最大，此时就不可能被分割，这就是为什么要令 K 为整个无向带权图中权值最大值，而其与顶点 \(T\) 相连的权值为 0 ，能量也为 0，此时一定会被分割。
\(B(A)\) 表示边界项，即上文中 n-links 边的权值，代表邻域像素点 \({p,q}\) 不连续的惩罚，当其差异越大时， \(B(A)\) 越趋近于 0，即当邻域像素点差异越大时，由上文 \(B_{\{p,q\}} \propto exp\left(-\frac{(I_p-I_q)^2}{2\sigma^2}\right) \cdot \frac{1}{dist(p,q)}\) 可知其权值越小，能量 \(B(A)\) 也越小，则越应该被分割。
参数 \(\lambda\) 用来表示区域项和边界项的重要性，初始计算时同样可以将其置为 1。
至于 \(\delta(A_p,A_q)\) 的个人理解为：当邻域像素点 \(p,q\) 被分配的二值变量不同时，即一个是前景另一个是背景，这时在进行能量计算时需要考虑连接亮点的边 \(\{p,q\}\) 的权值，否则可以不考虑其权值。

最后，本文割的具体实施是采用基于「最大流/最小割」算法进行的，而「最大流/最小割」算法针对的是有向图，所以需要把本文 S/T 图的每条无向边都转化为一去一回的两条有向边。

附录

　　该论文还有一个有意思的地方就是其证明了能量函数 E(A) 的值只与割 \(C\) 有关，最小割与最小化能量函数对应。证明过程如下：

由上文可知，对于图像中每个像素点，割 \(C\) 切断且仅切断一条 t-links 边，即要么切断 \(\{p,S\}\) 边，要么切断 \(\{p,T\}\) 边，毕竟一个像素点不可能同时属于前景和背景，也不可能既不属于前景也不属于背景；不属于同一终端顶点相连的 n-links 边也会被割 \(C\) 切断，因为如果不断开，则整个图还是相连的，而如果属于同一终端顶点相连的 n-links 边不应该被割 \(C\) 切断，否则，该割就不是最小割。
对于任意割 \(C\) ，可定义其对应的图像分割结果向量 \(A(C)\) ，如下： \[ A_p(C) = \begin{cases} "obj" & \text{if } \{p,T\} \in C \\ "bkg" & \text{if } \{p,S\} \in C \end{cases} \] 因此当 C 确定之后，分割结果 A 也就确定了。
利用上文给出的边的权值计算方法结合定义的 A 可得出割的代价（cost）： \[ \begin{align} |C| &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) + \sum_{p \in O} \lambda \cdot R_p(A_p(C)) + \sum_{p \in B} \lambda \cdot R_p(A_p(C)) + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \\ &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) +0 + 0 + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \end{align} \] 这里是因为 \(O \cap B = \emptyset\) ，而且如果确定像素点 p 属于前景或背景，则其割边的权值必为 0，相应的代价也为 0，eg：设 p 确定为前景，则必定断开 {p,T} 边，而 {p,T} 边的权值为 0。
又因为：\(E(A) = \lambda \cdot R(A) + B(A)\) ，则： \[ \begin{align} E(A(C)) &= \lambda \cdot R(A(C)) + B(A(C)) \\ &= \sum_{p \in P} \lambda \cdot R_p(A_p(C)) + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \\ &= \sum_{p \notin O \cup B} \lambda \cdot R_p(A_p(C)) + \sum_{p \in O} \lambda \cdot R_p("obj") + \sum_{p \in B} \lambda \cdot R_p("bkg") + \sum_{\{p.q\} \in N} B_{\{p,q\}} \cdot \delta((A_p(C)),A_q(C))) \end{align} \] 所以：\(|C| = E(A(C)) - \sum \limits_{p \in O} \lambda \cdot R_p("obj") - \sum \limits_{p \in B} \lambda \cdot R_p("bkg")\)
又因为对于标注确定的前景和背景，任意割 C 的 \(\sum \limits_{p \in O} \lambda \cdot R_p("obj")\) 和 \(\sum \limits_{p \in B} \lambda \cdot R_p("bkg")\) 都是完全确定不变的，可以令其为 const，则 \(E(A) = |C| + const\) ，当 \(C\) 取最小时，对应的 \(E(A)\) 也最小，即当割最小时，能量函数也最小。