----- -----《“GrabCut” — Interactive Foreground Extraction using Iterated Graph Cuts》

以下是我对grabcut的简单认识

grabcut是在graphcut基础上演变来的，grabcut算法利用了图像中的纹理（颜色）信息和边界（反差）信息，只要少量的用户交互操作即可得到比较好的分割结果。

一、Graph cut的介绍

首先，介绍一种图，如图1所示，这是一个特殊的图（graph），它的特殊之处在于除了中间的像素点之间（网格部分）构成图之外，他还有两个特殊顶点T（汇点）和S（源点），分别表示背景（Background）和前景（Object），以及各个像素和这两个点之间的连线，这两个特殊点是种子点，是用户在进行图像分割过程中标注出来的点（用户交互过程中指明哪里是背景，哪里是前景），如图2所示。

这种图的顶点确定了，那它各个点之间的边值又是怎么确定的呢？

这里我们就要引入最重要的公式了，利用这个公式不仅能解决上面的这个问题，还能帮我们进行图像分割，神奇，图像能量表示公式：

首先看公式（1）E(A)=lamda*R(A)+B(A),其中R(A)是区域项，表示像素A属于背景或者前景的概率，在图1中表示所有像素点与顶点T和顶点S的连线的边的权值，那我们怎么知道R(A) 为多少呢？这里要用到这个公式2，它等于该像素点属于前景或是背景的概率的负对数。

B(A)代表边界项，它对应图1中相邻像素点之间的边的权值，表示为两个像素点接近程度，这里用公式3来表示，两个像素点的像素值大小越相近，B（A）的值就越大，两个像素值大小越不相近，B(A)的值越小。

lamda表示R(A)和B(A)的比重，若lamda为0，代表只考虑B(A)，这里的lamda也是通过计算得出来的，具体详见第二部分。

这样我们不难发现，若是R(A)越小，代表像素属于前景或是背景的概率也就越大，B(A)越小，代表两个像素点不相近，不属于一类像素，应该分开，所以当E(A)能量越小，代表分割越准确。

那么我们如何进行分割呢？这里根据《Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images》这篇论文，提出了最小割方法。最小割算法原理可以参考https://blog.csdn.net/chinacoy/article/details/45040897。通俗易懂。

二、grabcut 介绍

与graph cut相比，grabcut具有以下优点：

1、只需一个长方形的目标框

2、增加少量的用户交互，分割更完美

3、border matting技术使目标分割边界更加自然

与graph cut 的不同之处在于：

1、Graph Cut 目标、背景是灰度直方图，Grab Cut是RGB三通道的混合高斯模型（GMM）

2、Graph Cut 分割一次完成，Grab Cut是不断进行分割估计和模型参数学习的迭代过程

3、Graph Cut 需要指定目标和背景种子点，Grab Cut 只需框选目标，允许不完全标注

上算法流程图，配合下面代码更容易看懂一些：

初始化部分模型迭代部分：

grabcut()函数内部代码：


/*M///////  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.////  By downloading, copying, installing or using the software you agree to this license.//  If you do not agree to this license, do not download, install,//  copy or use the software.//////                        Intel License Agreement//                For Open Source Computer Vision Library//// Copyright (C) 2000, Intel Corporation, all rights reserved.// Third party copyrights are property of their respective owners.//// Redistribution and use in source and binary forms, with or without modification,// are permitted provided that the following conditions are met:////   * Redistribution's of source code must retain the above copyright notice,//     this list of conditions and the following disclaimer.////   * Redistribution's in binary form must reproduce the above copyright notice,//     this list of conditions and the following disclaimer in the documentation//     and/or other materials provided with the distribution.////   * The name of Intel Corporation may not be used to endorse or promote products//     derived from this software without specific prior written permission.//// This software is provided by the copyright holders and contributors "as is" and// any express or implied warranties, including, but not limited to, the implied// warranties of merchantability and fitness for a particular purpose are disclaimed.// In no event shall the Intel Corporation or contributors be liable for any direct,// indirect, incidental, special, exemplary, or consequential damages// (including, but not limited to, procurement of substitute goods or services;// loss of use, data, or profits; or business interruption) however caused// and on any theory of liability, whether in contract, strict liability,// or tort (including negligence or otherwise) arising in any way out of// the use of this software, even if advised of the possibility of such damage.////M*/#include "precomp.hpp"#include "gcgraph.hpp"#include <limits>using namespace cv;/*This is implementation of image segmentation algorithm GrabCut described in"GrabCut — Interactive Foreground Extraction using Iterated Graph Cuts".Carsten Rother, Vladimir Kolmogorov, Andrew Blake.*//*GMM - Gaussian Mixture Model*/class GMM{public:static const int componentsCount = 5;GMM( Mat& _model );double operator()( const Vec3d color ) const;double operator()( int ci, const Vec3d color ) const;int whichComponent( const Vec3d color ) const;void initLearning();void addSample( int ci, const Vec3d color );void endLearning();private:void calcInverseCovAndDeterm( int ci );Mat model;double* coefs;double* mean;double* cov;double inverseCovs[componentsCount][3][3]; //协方差的逆矩阵double covDeterms[componentsCount];  //协方差的行列式double sums[componentsCount][3];double prods[componentsCount][3][3];int sampleCounts[componentsCount];int totalSampleCount;};//背景和前景各有一个对应的GMM（混合高斯模型）GMM::GMM( Mat& _model ){//一个像素的（唯一对应）高斯模型的参数个数或者说一个高斯模型的参数个数//一个像素RGB三个通道值，故3个均值，3*3个协方差，共用一个权值const int modelSize = 3/*mean*/ + 9/*covariance*/ + 1/*component weight*/;if( _model.empty() ){//一个GMM共有componentsCount个高斯模型，一个高斯模型有modelSize个模型参数_model.create( 1, modelSize*componentsCount, CV_64FC1 );_model.setTo(Scalar(0));}else if( (_model.type() != CV_64FC1) || (_model.rows != 1) || (_model.cols != modelSize*componentsCount) )CV_Error( CV_StsBadArg, "_model must have CV_64FC1 type, rows == 1 and cols == 13*componentsCount" );model = _model;//注意这些模型参数的存储方式：先排完componentsCount个coefs，再3*componentsCount个mean。//再3*3*componentsCount个cov。coefs = model.ptr<double>(0);  //GMM的每个像素的高斯模型的权值变量起始存储指针mean = coefs + componentsCount; //均值变量起始存储指针cov = mean + 3*componentsCount;  //协方差变量起始存储指针for( int ci = 0; ci < componentsCount; ci++ )if( coefs[ci] > 0 )//计算GMM中第ci个高斯模型的协方差的逆Inverse和行列式Determinant//为了后面计算每个像素属于该高斯模型的概率（也就是数据能量项）calcInverseCovAndDeterm( ci ); }//计算一个像素（由color=（B,G,R）三维double型向量来表示）属于这个GMM混合高斯模型的概率。//也就是把这个像素像素属于componentsCount个高斯模型的概率与对应的权值相乘再相加，//具体见论文的公式（10）。结果从res返回。//这个相当于计算Gibbs能量的第一个能量项（取负后）。double GMM::operator()( const Vec3d color ) const{double res = 0;for( int ci = 0; ci < componentsCount; ci++ )res += coefs[ci] * (*this)(ci, color );return res;}//计算一个像素（由color=（B,G,R）三维double型向量来表示）属于第ci个高斯模型的概率。//具体过程，即高阶的高斯密度模型计算式，具体见论文的公式（10）。结果从res返回double GMM::operator()( int ci, const Vec3d color ) const{double res = 0;if( coefs[ci] > 0 ){CV_Assert( covDeterms[ci] > std::numeric_limits<double>::epsilon() );Vec3d diff = color;double* m = mean + 3*ci;diff[0] -= m[0]; diff[1] -= m[1]; diff[2] -= m[2];double mult = diff[0]*(diff[0]*inverseCovs[ci][0][0] + diff[1]*inverseCovs[ci][1][0] + diff[2]*inverseCovs[ci][2][0])+ diff[1]*(diff[0]*inverseCovs[ci][0][1] + diff[1]*inverseCovs[ci][1][1] + diff[2]*inverseCovs[ci][2][1])+ diff[2]*(diff[0]*inverseCovs[ci][0][2] + diff[1]*inverseCovs[ci][1][2] + diff[2]*inverseCovs[ci][2][2]);res = 1.0f/sqrt(covDeterms[ci]) * exp(-0.5f*mult);}return res;}//返回这个像素最有可能属于GMM中的哪个高斯模型（概率最大的那个）int GMM::whichComponent( const Vec3d color ) const{int k = 0;double max = 0;for( int ci = 0; ci < componentsCount; ci++ ){double p = (*this)( ci, color );if( p > max ){k = ci;  //找到概率最大的那个，或者说计算结果最大的那个max = p;}}return k;}//GMM参数学习前的初始化，主要是对要求和的变量置零void GMM::initLearning(){for( int ci = 0; ci < componentsCount; ci++){sums[ci][0] = sums[ci][1] = sums[ci][2] = 0;prods[ci][0][0] = prods[ci][0][1] = prods[ci][0][2] = 0;prods[ci][1][0] = prods[ci][1][1] = prods[ci][1][2] = 0;prods[ci][2][0] = prods[ci][2][1] = prods[ci][2][2] = 0;sampleCounts[ci] = 0;}totalSampleCount = 0;}//增加样本，即为前景或者背景GMM的第ci个高斯模型的像素集（这个像素集是来用估//计计算这个高斯模型的参数的）增加样本像素。计算加入color这个像素后，像素集//中所有像素的RGB三个通道的和sums（用来计算均值），还有它的prods（用来计算协方差），//并且记录这个像素集的像素个数和总的像素个数（用来计算这个高斯模型的权值）。void GMM::addSample( int ci, const Vec3d color ){sums[ci][0] += color[0]; sums[ci][1] += color[1]; sums[ci][2] += color[2];prods[ci][0][0] += color[0]*color[0]; prods[ci][0][1] += color[0]*color[1]; prods[ci][0][2] += color[0]*color[2];prods[ci][1][0] += color[1]*color[0]; prods[ci][1][1] += color[1]*color[1]; prods[ci][1][2] += color[1]*color[2];prods[ci][2][0] += color[2]*color[0]; prods[ci][2][1] += color[2]*color[1]; prods[ci][2][2] += color[2]*color[2];sampleCounts[ci]++;totalSampleCount++;}//从图像数据中学习GMM的参数：每一个高斯分量的权值、均值和协方差矩阵；//这里相当于论文中“Iterative minimisation”的step 2void GMM::endLearning(){const double variance = 0.01;for( int ci = 0; ci < componentsCount; ci++ ){int n = sampleCounts[ci]; //第ci个高斯模型的样本像素个数if( n == 0 )coefs[ci] = 0;else{//计算第ci个高斯模型的权值系数coefs[ci] = (double)n/totalSampleCount; //计算第ci个高斯模型的均值double* m = mean + 3*ci;m[0] = sums[ci][0]/n; m[1] = sums[ci][1]/n; m[2] = sums[ci][2]/n;//计算第ci个高斯模型的协方差double* c = cov + 9*ci;c[0] = prods[ci][0][0]/n - m[0]*m[0]; c[1] = prods[ci][0][1]/n - m[0]*m[1]; c[2] = prods[ci][0][2]/n - m[0]*m[2];c[3] = prods[ci][1][0]/n - m[1]*m[0]; c[4] = prods[ci][1][1]/n - m[1]*m[1]; c[5] = prods[ci][1][2]/n - m[1]*m[2];c[6] = prods[ci][2][0]/n - m[2]*m[0]; c[7] = prods[ci][2][1]/n - m[2]*m[1]; c[8] = prods[ci][2][2]/n - m[2]*m[2];//计算第ci个高斯模型的协方差的行列式double dtrm = c[0]*(c[4]*c[8]-c[5]*c[7]) - c[1]*(c[3]*c[8]-c[5]*c[6]) + c[2]*(c[3]*c[7]-c[4]*c[6]);if( dtrm <= std::numeric_limits<double>::epsilon() ){//相当于如果行列式小于等于0，（对角线元素）增加白噪声，避免其变//为退化（降秩）协方差矩阵（不存在逆矩阵，但后面的计算需要计算逆矩阵）。// Adds the white noise to avoid singular covariance matrix.c[0] += variance;c[4] += variance;c[8] += variance;}//计算第ci个高斯模型的协方差的逆Inverse和行列式DeterminantcalcInverseCovAndDeterm(ci);}}}//计算协方差的逆Inverse和行列式Determinantvoid GMM::calcInverseCovAndDeterm( int ci ){if( coefs[ci] > 0 ){//取第ci个高斯模型的协方差的起始指针double *c = cov + 9*ci;double dtrm =covDeterms[ci] = c[0]*(c[4]*c[8]-c[5]*c[7]) - c[1]*(c[3]*c[8]-c[5]*c[6]) + c[2]*(c[3]*c[7]-c[4]*c[6]);//在C++中，每一种内置的数据类型都拥有不同的属性, 使用<limits>库可以获//得这些基本数据类型的数值属性。因为浮点算法的截断，所以使得，当a=2，//b=3时 10*a/b == 20/b不成立。那怎么办呢？//这个小正数（epsilon）常量就来了，小正数通常为可用给定数据类型的//大于1的最小值与1之差来表示。若dtrm结果不大于小正数，那么它几乎为零。//所以下式保证dtrm>0，即行列式的计算正确（协方差对称正定，故行列式大于0）。CV_Assert( dtrm > std::numeric_limits<double>::epsilon() );//三阶方阵的求逆inverseCovs[ci][0][0] =  (c[4]*c[8] - c[5]*c[7]) / dtrm;inverseCovs[ci][1][0] = -(c[3]*c[8] - c[5]*c[6]) / dtrm;inverseCovs[ci][2][0] =  (c[3]*c[7] - c[4]*c[6]) / dtrm;inverseCovs[ci][0][1] = -(c[1]*c[8] - c[2]*c[7]) / dtrm;inverseCovs[ci][1][1] =  (c[0]*c[8] - c[2]*c[6]) / dtrm;inverseCovs[ci][2][1] = -(c[0]*c[7] - c[1]*c[6]) / dtrm;inverseCovs[ci][0][2] =  (c[1]*c[5] - c[2]*c[4]) / dtrm;inverseCovs[ci][1][2] = -(c[0]*c[5] - c[2]*c[3]) / dtrm;inverseCovs[ci][2][2] =  (c[0]*c[4] - c[1]*c[3]) / dtrm;}}//计算beta，也就是Gibbs能量项中的第二项（平滑项）中的指数项的beta，用来调整//高或者低对比度时，两个邻域像素的差别的影响的，例如在低对比度时，两个邻域//像素的差别可能就会比较小，这时候需要乘以一个较大的beta来放大这个差别，//在高对比度时，则需要缩小本身就比较大的差别。//所以我们需要分析整幅图像的对比度来确定参数beta，具体的见论文公式（5）。/*Calculate beta - parameter of GrabCut algorithm.beta = 1/(2*avg(sqr(||color[i] - color[j]||)))*/static double calcBeta( const Mat& img ){double beta = 0;for( int y = 0; y < img.rows; y++ ){for( int x = 0; x < img.cols; x++ ){//计算四个方向邻域两像素的差别，也就是欧式距离或者说二阶范数//（当所有像素都算完后，就相当于计算八邻域的像素差了）Vec3d color = img.at<Vec3b>(y,x);if( x>0 ) // left  >0的判断是为了避免在图像边界的时候还计算，导致越界{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y,x-1);beta += diff.dot(diff);  //矩阵的点乘，也就是各个元素平方的和}if( y>0 && x>0 ) // upleft{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x-1);beta += diff.dot(diff);}if( y>0 ) // up{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x);beta += diff.dot(diff);}if( y>0 && x<img.cols-1) // upright{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x+1);beta += diff.dot(diff);}}}if( beta <= std::numeric_limits<double>::epsilon() )beta = 0;elsebeta = 1.f / (2 * beta/(4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2) ); //论文公式（5）return beta;}//计算图每个非端点顶点（也就是每个像素作为图的一个顶点，不包括源点s和汇点t）与邻域顶点//的边的权值。由于是无向图，我们计算的是八邻域，那么对于一个顶点，我们计算四个方向就行，//在其他的顶点计算的时候，会把剩余那四个方向的权值计算出来。这样整个图算完后，每个顶点//与八邻域的顶点的边的权值就都计算出来了。//这个相当于计算Gibbs能量的第二个能量项（平滑项），具体见论文中公式（4）/*Calculate weights of noterminal vertices of graph.beta and gamma - parameters of GrabCut algorithm.*/static void calcNWeights( const Mat& img, Mat& leftW, Mat& upleftW, Mat& upW, Mat& uprightW, double beta, double gamma ){//gammaDivSqrt2相当于公式（4）中的gamma * dis(i,j)^(-1)，那么可以知道，//当i和j是垂直或者水平关系时，dis(i,j)=1，当是对角关系时，dis(i,j)=sqrt(2.0f)。//具体计算时，看下面就明白了const double gammaDivSqrt2 = gamma / std::sqrt(2.0f);//每个方向的边的权值通过一个和图大小相等的Mat来保存leftW.create( img.rows, img.cols, CV_64FC1 );upleftW.create( img.rows, img.cols, CV_64FC1 );upW.create( img.rows, img.cols, CV_64FC1 );uprightW.create( img.rows, img.cols, CV_64FC1 );for( int y = 0; y < img.rows; y++ ){for( int x = 0; x < img.cols; x++ ){Vec3d color = img.at<Vec3b>(y,x);if( x-1>=0 ) // left  //避免图的边界{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y,x-1);leftW.at<double>(y,x) = gamma * exp(-beta*diff.dot(diff));}elseleftW.at<double>(y,x) = 0;if( x-1>=0 && y-1>=0 ) // upleft{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x-1);upleftW.at<double>(y,x) = gammaDivSqrt2 * exp(-beta*diff.dot(diff));}elseupleftW.at<double>(y,x) = 0;if( y-1>=0 ) // up{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x);upW.at<double>(y,x) = gamma * exp(-beta*diff.dot(diff));}elseupW.at<double>(y,x) = 0;if( x+1<img.cols && y-1>=0 ) // upright{Vec3d diff = color - (Vec3d)img.at<Vec3b>(y-1,x+1);uprightW.at<double>(y,x) = gammaDivSqrt2 * exp(-beta*diff.dot(diff));}elseuprightW.at<double>(y,x) = 0;}}}//检查mask的正确性。mask为通过用户交互或者程序设定的，它是和图像大小一样的单通道灰度图，//每个像素只能取GC_BGD or GC_FGD or GC_PR_BGD or GC_PR_FGD 四种枚举值，分别表示该像素//（用户或者程序指定）属于背景、前景、可能为背景或者可能为前景像素。具体的参考：//ICCV2001“Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images”//Yuri Y. Boykov Marie-Pierre Jolly /*Check size, type and element values of mask matrix.*/static void checkMask( const Mat& img, const Mat& mask ){if( mask.empty() )CV_Error( CV_StsBadArg, "mask is empty" );if( mask.type() != CV_8UC1 )CV_Error( CV_StsBadArg, "mask must have CV_8UC1 type" );if( mask.cols != img.cols || mask.rows != img.rows )CV_Error( CV_StsBadArg, "mask must have as many rows and cols as img" );for( int y = 0; y < mask.rows; y++ ){for( int x = 0; x < mask.cols; x++ ){uchar val = mask.at<uchar>(y,x);if( val!=GC_BGD && val!=GC_FGD && val!=GC_PR_BGD && val!=GC_PR_FGD )CV_Error( CV_StsBadArg, "mask element value must be equel""GC_BGD or GC_FGD or GC_PR_BGD or GC_PR_FGD" );}}}//通过用户框选目标rect来创建mask，rect外的全部作为背景，设置为GC_BGD，//rect内的设置为 GC_PR_FGD（可能为前景）/*Initialize mask using rectangular.*/static void initMaskWithRect( Mat& mask, Size imgSize, Rect rect ){mask.create( imgSize, CV_8UC1 );mask.setTo( GC_BGD );rect.x = max(0, rect.x);rect.y = max(0, rect.y);rect.width = min(rect.width, imgSize.width-rect.x);rect.height = min(rect.height, imgSize.height-rect.y);(mask(rect)).setTo( Scalar(GC_PR_FGD) );}//通过k-means算法来初始化背景GMM和前景GMM模型/*Initialize GMM background and foreground models using kmeans algorithm.*/static void initGMMs( const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM ){const int kMeansItCount = 10;  //迭代次数const int kMeansType = KMEANS_PP_CENTERS; //Use kmeans++ center initialization by Arthur and VassilvitskiiMat bgdLabels, fgdLabels; //记录背景和前景的像素样本集中每个像素对应GMM的哪个高斯模型，论文中的knvector<Vec3f> bgdSamples, fgdSamples; //背景和前景的像素样本集Point p;for( p.y = 0; p.y < img.rows; p.y++ ){for( p.x = 0; p.x < img.cols; p.x++ ){//mask中标记为GC_BGD和GC_PR_BGD的像素都作为背景的样本像素if( mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD )bgdSamples.push_back( (Vec3f)img.at<Vec3b>(p) );else // GC_FGD | GC_PR_FGDfgdSamples.push_back( (Vec3f)img.at<Vec3b>(p) );}}CV_Assert( !bgdSamples.empty() && !fgdSamples.empty() );//kmeans中参数_bgdSamples为：每行一个样本//kmeans的输出为bgdLabels，里面保存的是输入样本集中每一个样本对应的类标签（样本聚为componentsCount类后）Mat _bgdSamples( (int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0] );kmeans( _bgdSamples, GMM::componentsCount, bgdLabels,TermCriteria( CV_TERMCRIT_ITER, kMeansItCount, 0.0), 0, kMeansType );Mat _fgdSamples( (int)fgdSamples.size(), 3, CV_32FC1, &fgdSamples[0][0] );kmeans( _fgdSamples, GMM::componentsCount, fgdLabels,TermCriteria( CV_TERMCRIT_ITER, kMeansItCount, 0.0), 0, kMeansType );//经过上面的步骤后，每个像素所属的高斯模型就确定的了，那么就可以估计GMM中每个高斯模型的参数了。bgdGMM.initLearning();for( int i = 0; i < (int)bgdSamples.size(); i++ )bgdGMM.addSample( bgdLabels.at<int>(i,0), bgdSamples[i] );bgdGMM.endLearning();fgdGMM.initLearning();for( int i = 0; i < (int)fgdSamples.size(); i++ )fgdGMM.addSample( fgdLabels.at<int>(i,0), fgdSamples[i] );fgdGMM.endLearning();}//论文中：迭代最小化算法step 1：为每个像素分配GMM中所属的高斯模型，kn保存在Mat compIdxs中/*Assign GMMs components for each pixel.*/static void assignGMMsComponents( const Mat& img, const Mat& mask, const GMM& bgdGMM, const GMM& fgdGMM, Mat& compIdxs ){Point p;for( p.y = 0; p.y < img.rows; p.y++ ){for( p.x = 0; p.x < img.cols; p.x++ ){Vec3d color = img.at<Vec3b>(p);//通过mask来判断该像素属于背景像素还是前景像素，再判断它属于前景或者背景GMM中的哪个高斯分量compIdxs.at<int>(p) = mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD ?bgdGMM.whichComponent(color) : fgdGMM.whichComponent(color);}}}//论文中：迭代最小化算法step 2：从每个高斯模型的像素样本集中学习每个高斯模型的参数/*Learn GMMs parameters.*/static void learnGMMs( const Mat& img, const Mat& mask, const Mat& compIdxs, GMM& bgdGMM, GMM& fgdGMM ){bgdGMM.initLearning();fgdGMM.initLearning();Point p;for( int ci = 0; ci < GMM::componentsCount; ci++ ){for( p.y = 0; p.y < img.rows; p.y++ ){for( p.x = 0; p.x < img.cols; p.x++ ){if( compIdxs.at<int>(p) == ci ){if( mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD )bgdGMM.addSample( ci, img.at<Vec3b>(p) );elsefgdGMM.addSample( ci, img.at<Vec3b>(p) );}}}}bgdGMM.endLearning();fgdGMM.endLearning();}//通过计算得到的能量项构建图，图的顶点为像素点，图的边由两部分构成，//一类边是：每个顶点与Sink汇点t（代表背景）和源点Source（代表前景）连接的边，//这类边的权值通过Gibbs能量项的第一项能量项来表示。//另一类边是：每个顶点与其邻域顶点连接的边，这类边的权值通过Gibbs能量项的第二项能量项来表示。/*Construct GCGraph*/static void constructGCGraph( const Mat& img, const Mat& mask, const GMM& bgdGMM, const GMM& fgdGMM, double lambda,const Mat& leftW, const Mat& upleftW, const Mat& upW, const Mat& uprightW,GCGraph<double>& graph ){int vtxCount = img.cols*img.rows;  //顶点数，每一个像素是一个顶点int edgeCount = 2*(4*vtxCount - 3*(img.cols + img.rows) + 2);  //边数，需要考虑图边界的边的缺失//通过顶点数和边数创建图。这些类型声明和函数定义请参考gcgraph.hppgraph.create(vtxCount, edgeCount);Point p;for( p.y = 0; p.y < img.rows; p.y++ ){for( p.x = 0; p.x < img.cols; p.x++){// add nodeint vtxIdx = graph.addVtx();  //返回这个顶点在图中的索引Vec3b color = img.at<Vec3b>(p);// set t-weights           //计算每个顶点与Sink汇点t（代表背景）和源点Source（代表前景）连接的权值。//也即计算Gibbs能量（每一个像素点作为背景像素或者前景像素）的第一个能量项double fromSource, toSink;if( mask.at<uchar>(p) == GC_PR_BGD || mask.at<uchar>(p) == GC_PR_FGD ){//对每一个像素计算其作为背景像素或者前景像素的第一个能量项，作为分别与t和s点的连接权值fromSource = -log( bgdGMM(color) );toSink = -log( fgdGMM(color) );}else if( mask.at<uchar>(p) == GC_BGD ){//对于确定为背景的像素点，它与Source点（前景）的连接为0，与Sink点的连接为lambdafromSource = 0;toSink = lambda;}else // GC_FGD{fromSource = lambda;toSink = 0;}//设置该顶点vtxIdx分别与Source点和Sink点的连接权值graph.addTermWeights( vtxIdx, fromSource, toSink );// set n-weights  n-links//计算两个邻域顶点之间连接的权值。//也即计算Gibbs能量的第二个能量项（平滑项）if( p.x>0 ){double w = leftW.at<double>(p);graph.addEdges( vtxIdx, vtxIdx-1, w, w );}if( p.x>0 && p.y>0 ){double w = upleftW.at<double>(p);graph.addEdges( vtxIdx, vtxIdx-img.cols-1, w, w );}if( p.y>0 ){double w = upW.at<double>(p);graph.addEdges( vtxIdx, vtxIdx-img.cols, w, w );}if( p.x<img.cols-1 && p.y>0 ){double w = uprightW.at<double>(p);graph.addEdges( vtxIdx, vtxIdx-img.cols+1, w, w );}}}}//论文中：迭代最小化算法step 3：分割估计：最小割或者最大流算法/*Estimate segmentation using MaxFlow algorithm*/static void estimateSegmentation( GCGraph<double>& graph, Mat& mask ){//通过最大流算法确定图的最小割，也即完成图像的分割graph.maxFlow();Point p;for( p.y = 0; p.y < mask.rows; p.y++ ){for( p.x = 0; p.x < mask.cols; p.x++ ){//通过图分割的结果来更新mask，即最后的图像分割结果。注意的是，永远都//不会更新用户指定为背景或者前景的像素if( mask.at<uchar>(p) == GC_PR_BGD || mask.at<uchar>(p) == GC_PR_FGD ){if( graph.inSourceSegment( p.y*mask.cols+p.x /*vertex index*/ ) )mask.at<uchar>(p) = GC_PR_FGD;elsemask.at<uchar>(p) = GC_PR_BGD;}}}}//最后的成果：提供给外界使用的伟大的API：grabCut /*****参数说明：img——待分割的源图像，必须是8位3通道（CV_8UC3）图像，在处理的过程中不会被修改；mask——掩码图像，如果使用掩码进行初始化，那么mask保存初始化掩码信息；在执行分割的时候，也可以将用户交互所设定的前景与背景保存到mask中，然后再传入grabCut函数；在处理结束之后，mask中会保存结果。mask只能取以下四种值：GCD_BGD（=0），背景；GCD_FGD（=1），前景；GCD_PR_BGD（=2），可能的背景；GCD_PR_FGD（=3），可能的前景。如果没有手工标记GCD_BGD或者GCD_FGD，那么结果只会有GCD_PR_BGD或GCD_PR_FGD；rect——用于限定需要进行分割的图像范围，只有该矩形窗口内的图像部分才被处理；bgdModel——背景模型，如果为null，函数内部会自动创建一个bgdModel；bgdModel必须是单通道浮点型（CV_32FC1）图像，且行数只能为1，列数只能为13x5；fgdModel——前景模型，如果为null，函数内部会自动创建一个fgdModel；fgdModel必须是单通道浮点型（CV_32FC1）图像，且行数只能为1，列数只能为13x5；iterCount——迭代次数，必须大于0；mode——用于指示grabCut函数进行什么操作，可选的值有：GC_INIT_WITH_RECT（=0），用矩形窗初始化GrabCut；GC_INIT_WITH_MASK（=1），用掩码图像初始化GrabCut；GC_EVAL（=2），执行分割。*/void cv::grabCut( InputArray _img, InputOutputArray _mask, Rect rect,InputOutputArray _bgdModel, InputOutputArray _fgdModel,int iterCount, int mode ){Mat img = _img.getMat();Mat& mask = _mask.getMatRef();Mat& bgdModel = _bgdModel.getMatRef();Mat& fgdModel = _fgdModel.getMatRef();if( img.empty() )CV_Error( CV_StsBadArg, "image is empty" );if( img.type() != CV_8UC3 )CV_Error( CV_StsBadArg, "image mush have CV_8UC3 type" );GMM bgdGMM( bgdModel ), fgdGMM( fgdModel );Mat compIdxs( img.size(), CV_32SC1 );if( mode == GC_INIT_WITH_RECT || mode == GC_INIT_WITH_MASK ){if( mode == GC_INIT_WITH_RECT )initMaskWithRect( mask, img.size(), rect );else // flag == GC_INIT_WITH_MASKcheckMask( img, mask );initGMMs( img, mask, bgdGMM, fgdGMM );}if( iterCount <= 0)return;if( mode == GC_EVAL )checkMask( img, mask );const double gamma = 50;const double lambda = 9*gamma;const double beta = calcBeta( img );Mat leftW, upleftW, upW, uprightW;calcNWeights( img, leftW, upleftW, upW, uprightW, beta, gamma );for( int i = 0; i < iterCount; i++ ){GCGraph<double> graph;assignGMMsComponents( img, mask, bgdGMM, fgdGMM, compIdxs );learnGMMs( img, mask, compIdxs, bgdGMM, fgdGMM );constructGCGraph(img, mask, bgdGMM, fgdGMM, lambda, leftW, upleftW, upW, uprightW, graph );estimateSegmentation( graph, mask );}}

添加个人理解部分：

一、类型转化

InputArray、InputOutputArray、InputOutputArray、InputOutputArray转化为Mat类型

二、检查输入图像的合法性

是否为空、是否是CV_8UC3类型

三、建立前景、背景高斯混合模型（GMM）

确定每个高斯模型参数个数=均值（3）+协方差（3*3）+权值（1）

建立混合高斯模型（参数=模型个数*每个模型的参数）

定义模型参数存储方式：权值、均值、协方差变量的起始指针

四、构造Mat 类型的单通道compIdxs（与img大小相同）,用于存储每个像素属于哪个高斯模型

五、判断grabcut函数进行的操作

GC_INIT_WITH_RECT：（初始矩形框）创建一个cv_8uc1类型的和图片大小一样的mask,设置为背景；在这个mask上选择和矩形一样大小地方设置为前景

GC_INIT_WITH_MASK：（检查mask合法性）mask是否为空、是否是cv_8uc1类型，行列是否与图片大小一致、mask取值是否为4中枚举值

六、初始化GMMS：

确定k-means迭代次数(10),类型(KMEANS_PP_CENTERS).

扫描像素点，根据mask将像素点加入前景、背景样本集

执行k-means聚类，输出每一个样本（像素）对应哪个高斯模型

bdgGMM. initLearning()：参数、求和变量置零

bdgGMM. addSample( int ci, const Vec3d color )：计算每个高斯模型color[i](i=1,2,3)的总和，协方差，第ci高斯模型样本个数，混合高斯模型总样本个数

bdgGMM. endLearning():计算第ci高斯模型的权值、均值、协方差、协方差的逆和行列式

七、计算每个边的权重：

初始化gama=50,lamda=9*gama(无穷大)，计算bata(欧氏距离)

利用公式gamma * exp(-beta*diff.dot(diff))计算像素四个方向的权重

权重保存在mat理性的变量

八、fdgGMM同理

九、迭代过程

step1:为每个像素分配模型

1、通过mask判断像素属于前景还是背景

2、判断属于哪个高斯分量（高斯模型密度计算式）

3、将概率最大的ci存储到对应的compidx中

step2:学习高斯模型参数

1、利用 bdgGMM. initLearning()： bdgGMM. addSample( int ci, const Vec3d color ),bdgGMM. endLearning(),计算参数。

2、构建图：每个点与源点和汇点的连接权值、邻域顶点间的连接权值

step3:分割

1、最大流算法分割成两部分

2、更新mask（不更新指定的前景和背景）

回到step1直到迭代结束。

应用grabcut进行图片分割：

#pragma once
#include"pch.h"
#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"
//#include "BorderMatting.h"#include <iostream>using namespace std;using namespace cv;static void help(){cout << "\nThis program demonstrates GrabCut segmentation -- select an object in a region\n""and then grabcut will attempt to segment it out.\n""Call:\n""./grabcut <image_name>\n""\nSelect a rectangular area around the object you want to segment\n" <<"\nHot keys: \n""\tESC - quit the program\n""\tr - restore the original image\n""\tn - next iteration\n""\n""\tleft mouse button - set rectangle\n""\n""\tCTRL+left mouse button - set GC_BGD pixels\n""\tSHIFT+left mouse button - set CG_FGD pixels\n""\n""\tCTRL+right mouse button - set GC_PR_BGD pixels\n""\tSHIFT+right mouse button - set CG_PR_FGD pixels\n" << endl;}const Scalar RED = Scalar(0, 0, 255);const Scalar PINK = Scalar(230, 130, 255);const Scalar BLUE = Scalar(255, 0, 0);const Scalar LIGHTBLUE = Scalar(255, 255, 160);const Scalar GREEN = Scalar(0, 255, 0);const int BGD_KEY = CV_EVENT_FLAG_CTRLKEY;  //Ctrl键const int FGD_KEY = CV_EVENT_FLAG_SHIFTKEY; //Shift键static void getBinMask(const Mat& comMask, Mat& binMask){if (comMask.empty() || comMask.type() != CV_8UC1)CV_Error(CV_StsBadArg, "comMask is empty or has incorrect type (not CV_8UC1)");if (binMask.empty() || binMask.rows != comMask.rows || binMask.cols != comMask.cols)binMask.create(comMask.size(), CV_8UC1);binMask = comMask & 1;  //得到mask的最低位,实际上是只保留确定的或者有可能的前景点当做mask}class GCApplication{public:enum { NOT_SET = 0, IN_PROCESS = 1, SET = 2 };static const int radius = 2;static const int thickness = -1;void reset();void setImageAndWinName(const Mat& _image, const string& _winName);void showImage() const;void mouseClick(int event, int x, int y, int flags, void* param);int nextIter();int getIterCount() const { return iterCount; }private:void setRectInMask();void setLblsInMask(int flags, Point p, bool isPr);const string* winName;const Mat* image;Mat mask;Mat bgdModel, fgdModel;uchar rectState, lblsState, prLblsState;bool isInitialized;Rect rect;vector<Point> fgdPxls, bgdPxls, prFgdPxls, prBgdPxls;int iterCount;
//  BorderMatting };/*给类的变量赋值*/void GCApplication::reset(){if (!mask.empty())mask.setTo(Scalar::all(GC_BGD));bgdPxls.clear(); fgdPxls.clear();prBgdPxls.clear();  prFgdPxls.clear();isInitialized = false;rectState = NOT_SET;    //NOT_SET == 0lblsState = NOT_SET;prLblsState = NOT_SET;iterCount = 0;}/*给类的成员变量赋值而已*/void GCApplication::setImageAndWinName(const Mat& _image, const string& _winName){if (_image.empty() || _winName.empty())return;image = &_image;winName = &_winName;mask.create(image->size(), CV_8UC1);reset();}/*显示4个点，一个矩形和图像内容，因为后面的步骤很多地方都要用到这个函数，所以单独拿出来*/void GCApplication::showImage() const{if (image->empty() || winName->empty())return;Mat res;Mat binMask;if (!isInitialized)image->copyTo(res);else{getBinMask(mask, binMask);image->copyTo(res, binMask);  //按照最低位是0还是1来复制，只保留跟前景有关的图像，比如说可能的前景，可能的背景}vector<Point>::const_iterator it;/*下面4句代码是将选中的4个点用不同的颜色显示出来*/for (it = bgdPxls.begin(); it != bgdPxls.end(); ++it)  //迭代器可以看成是一个指针circle(res, *it, radius, BLUE, thickness);for (it = fgdPxls.begin(); it != fgdPxls.end(); ++it)  //确定的前景用红色表示circle(res, *it, radius, RED, thickness);for (it = prBgdPxls.begin(); it != prBgdPxls.end(); ++it)circle(res, *it, radius, LIGHTBLUE, thickness);for (it = prFgdPxls.begin(); it != prFgdPxls.end(); ++it)circle(res, *it, radius, PINK, thickness);/*画矩形*/if (rectState == IN_PROCESS || rectState == SET)rectangle(res, Point(rect.x, rect.y), Point(rect.x + rect.width, rect.y + rect.height), GREEN, 2);imshow(*winName, res);}/*该步骤完成后，mask图像中rect内部是3，外面全是0*/void GCApplication::setRectInMask(){assert(!mask.empty());mask.setTo(GC_BGD);   //GC_BGD == 0rect.x = max(0, rect.x);rect.y = max(0, rect.y);rect.width = min(rect.width, image->cols - rect.x);rect.height = min(rect.height, image->rows - rect.y);(mask(rect)).setTo(Scalar(GC_PR_FGD));    //GC_PR_FGD == 3，矩形内部,为可能的前景点}void GCApplication::setLblsInMask(int flags, Point p, bool isPr){vector<Point> *bpxls, *fpxls;uchar bvalue, fvalue;if (!isPr) //确定的点{bpxls = &bgdPxls;fpxls = &fgdPxls;bvalue = GC_BGD;    //0fvalue = GC_FGD;    //1}else    //概率点{bpxls = &prBgdPxls;fpxls = &prFgdPxls;bvalue = GC_PR_BGD; //2fvalue = GC_PR_FGD; //3}if (flags & BGD_KEY){bpxls->push_back(p);circle(mask, p, radius, bvalue, thickness);   //该点处为2}if (flags & FGD_KEY){fpxls->push_back(p);circle(mask, p, radius, fvalue, thickness);   //该点处为3}}/*鼠标响应函数，参数flags为CV_EVENT_FLAG的组合*/void GCApplication::mouseClick(int event, int x, int y, int flags, void*){// TODO add bad args checkswitch (event){case CV_EVENT_LBUTTONDOWN: // set rect or GC_BGD(GC_FGD) labels左键按下{bool isb = (flags & BGD_KEY) != 0,isf = (flags & FGD_KEY) != 0;if (rectState == NOT_SET && !isb && !isf)//只有左键按下时{rectState = IN_PROCESS; //表示正在画矩形rect = Rect(x, y, 1, 1);}if ((isb || isf) && rectState == SET) //按下了alt键或者shift键，且画好了矩形，表示正在画前景背景点lblsState = IN_PROCESS;}break;case CV_EVENT_RBUTTONDOWN: // set GC_PR_BGD(GC_PR_FGD) labels右键按下{bool isb = (flags & BGD_KEY) != 0,isf = (flags & FGD_KEY) != 0;if ((isb || isf) && rectState == SET) //正在画可能的前景背景点prLblsState = IN_PROCESS;}break;case CV_EVENT_LBUTTONUP://左键松起if (rectState == IN_PROCESS){rect = Rect(Point(rect.x, rect.y), Point(x, y));   //矩形结束rectState = SET;setRectInMask();assert(bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty());showImage();}if (lblsState == IN_PROCESS)   //已画了前后景点{setLblsInMask(flags, Point(x, y), false);    //画出前景点lblsState = SET;showImage();}break;case CV_EVENT_RBUTTONUP://右键松起if (prLblsState == IN_PROCESS){setLblsInMask(flags, Point(x, y), true); //画出背景点prLblsState = SET;showImage();}break;case CV_EVENT_MOUSEMOVE:if (rectState == IN_PROCESS){rect = Rect(Point(rect.x, rect.y), Point(x, y));assert(bgdPxls.empty() && fgdPxls.empty() && prBgdPxls.empty() && prFgdPxls.empty());showImage();    //不断的显示图片}else if (lblsState == IN_PROCESS){setLblsInMask(flags, Point(x, y), false);showImage();}else if (prLblsState == IN_PROCESS){setLblsInMask(flags, Point(x, y), true);showImage();}break;}}/*该函数进行grabcut算法，并且返回算法运行迭代的次数*/int GCApplication::nextIter(){if (isInitialized)//使用grab算法进行一次迭代，参数2为mask，里面存的mask位是：矩形内部除掉那些可能是背景或者已经确定是背景后的所有的点，且mask同时也为输出//保存的是分割后的前景图像grabCut(*image, mask, rect, bgdModel, fgdModel, 1);else{if (rectState != SET)return iterCount;if (lblsState == SET || prLblsState == SET)grabCut(*image, mask, rect, bgdModel, fgdModel, 1, GC_INIT_WITH_MASK);elsegrabCut(*image, mask, rect, bgdModel, fgdModel, 1, GC_INIT_WITH_RECT);isInitialized = true;}iterCount++;bgdPxls.clear(); fgdPxls.clear();prBgdPxls.clear(); prFgdPxls.clear();return iterCount;}GCApplication gcapp;static void on_mouse(int event, int x, int y, int flags, void* param){gcapp.mouseClick(event, x, y, flags, param);}int main(int argc, char** argv){string filename = "D:/picture/source/1.jpg";Mat image = imread(filename, 1);if (image.empty()){cout << "\n Durn, couldn't read image filename " << filename << endl;return 1;}help();const string winName = "image";cvNamedWindow(winName.c_str(), CV_WINDOW_AUTOSIZE);cvSetMouseCallback(winName.c_str(), on_mouse, 0);//window_name 回掉函数需要注册到的窗口，即产生事件的窗口,on_mouse 在注册窗口点击鼠标时，执行的回掉函数,param 用于传递到回掉函数的参数gcapp.setImageAndWinName(image, winName);gcapp.showImage();for (;;){int c = cvWaitKey(0);switch ((char)c){case '\x1b'://esc退出cout << "Exiting ..." << endl;goto exit_main;case 'r':cout << endl;gcapp.reset();gcapp.showImage();break;case 'n':int iterCount = gcapp.getIterCount();cout << "<" << iterCount << "... ";int newIterCount = gcapp.nextIter();if (newIterCount > iterCount){gcapp.showImage();cout << iterCount << ">" << endl;}elsecout << "rect must be determined>" << endl;break;}}exit_main:cvDestroyWindow(winName.c_str());return 0;}

分割结果显示：

参考https://blog.csdn.net/zouxy09/article/details/8535087

https://blog.csdn.net/wh1319501722/article/details/80272167

https://blog.csdn.net/yangyangyang20092010/article/details/25967303

迭代图像切割技术的交互式前景提取相关推荐

图像分割与提取：交互式前景提取（附OpenCV代码实现）
一.简介经典的前景提取技术主要使用纹理(颜色)信息,如魔术棒工具,或根据边缘(对比度)信息,如智能剪刀等完成.2004 年,微软研究院(剑桥)的 Rother 等人在论文 GrabCut: Inte ...
OpenCV系列之交互式前景提取使用GrabCut算法 | 三十五
目标在本章中, 我们将看到GrabCut算法来提取图像中的前景我们将为此创建一个交互式应用程序. 理论 GrabCut算法由英国微软研究院的Carsten Rother,Vladimir Kolm ...
opencv29:使用Grabcut算法的交互式前景提取
目标在这一章当中看到 GrabCut算法来提取图像中的前景为此创建一个交互式应用程序理论 GrabCut 算法由英国剑桥微软研究院的 Carsten Rother.Vladimir Kolmo ...
OpenCV使用 GrabCut 算法进行交互式前景提取
OpenCV使用 GrabCut 算法进行交互式前景提取 1. 效果图 2. 源码参考这篇博客将介绍如何使用Python,OpenCV中的GrabCut 算法来提取图像中的前景,并为此创建一个交互 ...
OpenCV中的图像处理 —— 霍夫线 / 圈变换 + 图像分割（分水岭算法） + 交互式前景提取（GrabCut算法）
OpenCV中的图像处理 -- 霍夫线 / 圈变换 + 图像分割(分水岭算法) + 交互式前景提取(GrabCut算法)
【OpenCV-Python】教程：3-16 利用Grabcut交互式前景提取
OpenCV Python Grabcut分割 [目标] Grabcut 算法创建一个交互程序 [理论] 从用户角度是如何工作的呢?用户在需要的目标上初始绘制一个矩形,前景目标必须完全在矩形内部,算 ...
Python OpenCV 交互式前景提取自动抠图
这是需要抠图的原图像,文件名为 "messi5.jpg" 使用矩形方式处理: # -*- coding: utf-8 -*-import numpy as np import cv ...
opencv入门：使用交互式进行前景提取
交互式前景提取在开始提取前景之前,先用一个矩形框指定前景区域所在的大致位置范围,然后不断的迭代的分割,直到达到最好的效果,经过上述处理,结果可能不理想,存在前景没有提取,或背景被提取,此时需要用户干 ...
［OpenCV_GrubCut］实现交互式图像分割提取前景--Python抠图
这部分内容是几个月前做的项目,一直没时间整理记录,在这里随便写一下方便日后回忆．＂GrabCut＂:使用迭代图形切割的交互式前景提取工具,用于在分割任务中按像素标记图像数据． OpenCV官网例子 ...

迭代图像切割技术的交互式前景提取

一、Graph cut的介绍

二、grabcut 介绍

迭代图像切割技术的交互式前景提取相关推荐

最新文章

热门文章