经典图割算法中图的构建及实现:one-Cut

经典图割算法中图的构建及实现之one-cut-1

本文目的：

讲解目前典型的3种图割算法：graph-cut、grab-but、one-cut。本文主要讲解one-cut的方法在应用时，准则函数与图构建关系，如何构建图，以及如何代码实现图的构建。图割的原理网上文章和论文已介绍比较详细，不再详细介绍。One-cut的人工干预图像分割算法有2种形式：人工设置正负种子点(graph-cut的方式)；人工设置背景区(grab-cut画框的方式)

• one-cut思想

one-cut最主要的改进在于将原始图割的先验惩罚项进行了修改。使用全图的前景与背景直方图的L1距离表示先验惩罚项。作者表示该修改能提升性能，并且解决了grab-cut的画框方式中NP难的问题。在grab-cut中需要用EM算法迭代，其不一定能求解到全局最优解，而one-cut可以直接求解到准则函数的最优解。该篇文章介绍选择正负种子点的one-cut算法，设置背景区方式(bounding-box)的one-cut在下一篇介绍。

• 准则函数

                                     

准则函数如上，相比graph-cut，先验项有所不同。此处直接使用直方图的L1距离为惩罚项。作者的论文中有相关相应的对比实验，说明该方法的优势。

       我们可以看到有如下关系成立：

                                                      

       K表示直方图bin数量，表示正负类各自在直方图第K个bin中的数量，为全集。

1.术语：

 

1. 与S和T链接的边叫t-link(红线与绿线)，领域之间的链接边叫n-link(黑线)，与辅节点相连的a-link(黄线)。其中红线进一步称为s-t-link，绿线进一步称为t-t-link，仅有种子点(红色为正种子点、绿色为负种子点)才存在t-link。
2. 黄色点表示的是辅节点，该节点的用处是构造R(A)项，有多少个直方图Bin就有多少个辅节点。属于同一bin的像素点连接向同一个辅节点。
3. 黑线的权值对应的是B(A)项，黄线的权值对应的是R(A)项，而红线和绿线对应的种子点硬约束项。
4. 权值用w表示。
5. 蓝色节点表示类别标志节点，S表示正类类标节点，T表示负类类标节点。

       最终通过求最小割之后，与节点S相连的所有红色、白色节点(图像像素点)属于一类，同理与节点T相连的所有绿色、白色节点属于另一类。两类被最小割割开，割值即是准则函数值。

2.图的建立

       拿到待分割的图像后，图的节点与边已确定，即图的形状已确定下来。仅仅需要做的就是给图中所有边赋相应的权值，如下图所示。

       图中的边有3种情况：种子点的t-link；每个像素点对应的a-link；像素领域关系的n-link。接下来将说明每一种边的权值取值。

1).种子点t-link权值：种子点认为是硬约束，其用户预设2类种子点后，种子点类别不会随分割算法而改变。

a.对于正类别种子点，s-t-link必须保留，t-t-link必须割去。工程中，通过将s-t-link权值设置为超级大值，t-t-link设置为0。保证一定仅仅割去t-t-link，否则一定不是最小割，因为当前w(s-t-link)权值是超级大值，割去这条边的代价一定是最大的。

b.反之同理。

2).a-link权值：原始图像确定后其颜色直方图即确定。将所有像素点按照其颜色值所属建立与对应颜色Bin节点的连接。为撒该节点的权值对应R(A)项，请各位读者自己任意取一个分割看看，其中常数项 直接去掉，因为其不影响最值。

3).n-link权值：n-link用于度量相邻像素点之间颜色的差异性。设一对相邻点Pi，Pj，则n-link(Pi-Pj)的权值w等于：

                                                                 

         其中，dist()是距离函数，表示点之间的图像距离。即4领域下，所以领域点距离均为1；8领域下，对角像素点距离为 ；在5*5领域下，对角像素点距离为 。

         设种子点的超级大值是1000， =1，K=4(灰度空间0-255被分割成4个Bin)， ， 。图像是3*2的灰度图，数字表示灰度值，红色和蓝色节点表示用户选择的正负种子点。

       如上图所示，将所有边的权值赋值后，图就建立完毕。剩余则直接运用最小割算法即可求解。最小割算法有很多，包括graph-cut作者提出的快速算法An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision。Opencv即采用该算法计算最小割。

3.建立图的代码实现：

1)图的构建：

    for (int i = 0; i<inputImg.rows; i++)
	{
		for (int j = 0; j<inputImg.cols; j++)
		{
			// this is the node id for the current pixel
			GraphType::node_id currNodeId = i * inputImg.cols + j;

			// add hard constraints based on scribbles
			if (fgScribbleMask.at<uchar>(i, j) == 255)
				myGraph->add_tweights(currNodeId, (int)ceil(INT32_CONST * HARD_CONSTRAINT_CONST + 0.5), 0);
			else if (bgScribbleMask.at<uchar>(i, j) == 255)
				myGraph->add_tweights(currNodeId, 0, (int)ceil(INT32_CONST * HARD_CONSTRAINT_CONST + 0.5));

			// You can now access the pixel value with cv::Vec3b
			float b = (float)inputImg.at<Vec3b>(i, j)[0];
			float g = (float)inputImg.at<Vec3b>(i, j)[1];
			float r = (float)inputImg.at<Vec3b>(i, j)[2];

			// go over the neighbors
			for (int si = -NEIGHBORHOOD; si <= NEIGHBORHOOD; si++)
			{
				int ni = i + si;
				// outside the border - skip
				if (ni < 0 || ni >= inputImg.rows)
					continue;

				for (int sj = 0; sj <= NEIGHBORHOOD; sj++)
				{
					int nj = j + sj;
					// outside the border - skip
					if (nj < 0 || nj >= inputImg.cols)
						continue;

					// same pixel - skip
					// down pointed edge, this edge will be counted as an up edge for the other pixel
					if (si >= 0 && sj == 0)
						continue;

					// diagonal exceed the radius - skip
					if ((si*si + sj*sj) > NEIGHBORHOOD*NEIGHBORHOOD)
						continue;


					// this is the node id for the neighbor
					GraphType::node_id nNodeId = (i + si) * inputImg.cols + (j + sj);

					float nb = (float)inputImg.at<Vec3b>(i + si, j + sj)[0];
					float ng = (float)inputImg.at<Vec3b>(i + si, j + sj)[1];
					float nr = (float)inputImg.at<Vec3b>(i + si, j + sj)[2];

					//   ||I_p - I_q||^2  /   2 * sigma^2
					float currEdgeStrength = exp(-((b - nb)*(b - nb) + (g - ng)*(g - ng) + (r - nr)*(r - nr)) / (2 * varianceSquared));
					//float currEdgeStrength = 0;
					float currDist = sqrt((float)si*(float)si + (float)sj*(float)sj);

					// this is the edge between the current two pixels (i,j) and (i+si, j+sj)
					currEdgeStrength = ((float)EDGE_STRENGTH_WEIGHT * currEdgeStrength + (float)(1 - EDGE_STRENGTH_WEIGHT)) / currDist;
					int edgeCapacity = /* capacities */ (int)ceil(INT32_CONST*currEdgeStrength + 0.5);
					//edgeCapacity = 0;
					myGraph->add_edge(currNodeId, nNodeId, edgeCapacity, edgeCapacity);

				}
			}
			// add the adge to the auxiliary node
			int currBin = (int)binPerPixelImg.at<float>(i, j);

			myGraph->add_edge(currNodeId, (GraphType::node_id)(currBin + inputImg.rows * inputImg.cols),
				/* capacities */ (int)ceil(INT32_CONST*bha_slope + 0.5), (int)ceil(INT32_CONST*bha_slope + 0.5));
		}//for (int j = 0; j<inputImg.cols; j++)

	}//for (int i = 0; i<inputImg.rows; i++)

2)分割效果：

                                                        

       从实验结果来看，图中如果有多个类，该方法一般不能取得较好结果。对于2类的图像，该方法效果很好，最后仅需再加上一些空洞填补、小区域过滤等操作就好。

3）完整资源代码：https://download.csdn.net/download/jy02660221/10776581，最后鄙视一下CSDN最低下载分数限制。