grab cut算法是graph cut算法的改进。在理解grab cut算之前，应该学习一下graph cut算法的概念及实现方式。

我搜集了一些graph cut资料：http://yunpan.cn/QGDVdBXwkXutH

grab cut算法详细描述见资料中的pdf文件：“GrabCut” — Interactive Foreground Extraction using Iterated Graph Cuts

grab cut算法是一种基于图论的图像分割方法，首先要定义一个Gibbs能量函数，然后求解这个函数的min-cut，这个min-cut就是前景背景的分割像素集合。

1. 能量函数的定义

在grab cut算法中，能量函数定义为：

其中U函数部分表示能量函数的区域数据项，V函数表示能量函数的光滑项(边界项)。

我们使用混合多高斯模型D(x)表示某个像素属于前景或背景的概率，这里K=5，表示第i个单高斯函数对概率贡献的权重系数，所以有， 为第i单高斯函数。为第i个单高斯函数的均值，为第i个单高斯函数的协方差。

单高斯函数g公式为：

区域数据项U函数为：

这里n表示图像中的第n个像素。【注：这儿的，是对前面的D(x)取负对数】：

grab cut算法的输入图像是RGB 3通道的图像,对于输入图像，我们用两个混合多高斯模型来分别表示前景和背景。

光滑项V函数为：

C是相邻颜色对的集合，是一个常量值50，

这些公式直接理解有点困难，下面我们结合程序看看grabcut算法中，如何计算能量公式，以及如何分段。

1、首先定义混合多高斯模型

class GMM
{
public:
    static const int componentsCount = 5; //对应K=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]; //所有样本bgr三个颜色分量的和，用来计算权重系数
    double prods[componentsCount][3][3]; //所有样本bgr颜色的行列式值，用来计算协方差
    int sampleCounts[componentsCount]; //每个单高斯函数的样本数
    int totalSampleCount;  //样本总数

};

从上面混合多高斯公式可以知道，只要确定了三个参数：权重系数、均值、协方差，就可以根据当前像素点的bgr值确定当前像素属于前景和背景的概率D(x)，所以在GMM类中，我们定义三个指针，分别表示权重系数，均值和协方差。因为当前像素用bgr值表示，所以均值其实为3个double数，再加上K=5(5个单高斯函数组成的多高斯混合函数)，总共15双精度值，而权重系数则为5个双精度值，cov公共3*3*5=45个双精度值。

double* coefs; //权重系数
double* mean; //均值
double* cov; //协方差

在GMM的构造函数中，我们会创建一个1维的矩阵，总共65个双精度数，权重系数指向矩阵数据头，均值指向第6个双精度数，协方差指向第21个双精度数。

_model.create( 1, modelSize*componentsCount, CV_64FC1 );
    _model.setTo(Scalar(0));

coefs = model.ptr<double>(0);
mean = coefs + componentsCount;
cov = mean + 3*componentsCount;

2.初始化GMM变量

我们定义了两个变量

GMM bgdGMM( bgdModel ), fgdGMM( fgdModel );

分别表示和前景的混合多高斯模型。

首先我们根据选定的四边形框来初始化mask图像，四边形框外的像素是背景，值为GC_BGD ，四边形内的像素可能是前景，值为GC_PR_FGD。

/*
  Initialize mask using rectangular.
  设置mask的初始值，四边形框内圈定的像素值为GC_PR_FGD
*/
void gGrabCut::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) );
}

之后，我们会根据mask图像，读入样本数据。前景GMM的样本数据放在变量fgdSamples中，背景GMM的样本数据放入变量bgdSamples中。fgdSamples和bgdSamples中存放得都是一些bgr颜色值。

Mat bgdLabels, fgdLabels;
vector<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++ )
    {
        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_FGD
            fgdSamples.push_back( (Vec3f)img.at<Vec3b>(p) );
    }
}
     之后我们会根据kmeans聚类算法，计算得到当前像素属于前景或背景混合多高斯变量中的第几个单高斯函数，结果放在bgdSamples, fgdSamples中，值为0-4。

//一行是一个数据样本，3列是b，g，r三个属性
    Mat _bgdSamples( (int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0] );
    //GMM::componentsCount聚类的个数
    //KMEANS_PP_CENTERS是采用Arthur & Vassilvitskii (2007) k-means++: The Advantages of Careful Seeding获取初始化种子点
    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 );

下面我们计算得到均值、协方差和权重系数：

权重系数为属于某个单高斯函数的采样像素数量除以所有采样像素的数量。每个单高斯函数的均值为所有属于该函数的采样像素颜色和除以属于该函数的颜色采样数量。协方差的公式比较复杂，大家看看下面代码中c[0]-c[9]的计算就ok了。注意其中的函数calcInverseCovAndDeterm用来计算协方差矩阵的行列式值以及逆矩阵，这些值在计算能量公式的数据项函数U时候使用。

//计算得到均值、协方差以及权重系数
void GMM::endLearning()
{
    const double variance = 0.01;
    for( int ci = 0; ci < componentsCount; ci++ )
    {
        int n = sampleCounts[ci];
        if( n == 0 )
            coefs[ci] = 0;
        else
        {
            coefs[ci] = (double)n/totalSampleCount;

double* m = mean + 3*ci;
            m[0] = sums[ci][0]/n; m[1] = sums[ci][1]/n; m[2] = sums[ci][2]/n;

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];

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() )
            {
               // Adds the white noise to avoid singular covariance matrix.
                c[0] += variance;
                c[4] += variance;
                c[8] += variance;
            }

calcInverseCovAndDeterm(ci);
        }
    }
}

下面我们开始计算能量公式中光滑性函数V：

（注：=1)，

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 );

我们会在beta函数计算公式中的beta值，其中4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2为邻接距离的数量。

/*
  计算光滑性函数中的beta值
  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
            {
                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;
    else //除以邻接距离的数量
        beta = 1.f / (2 * beta/(4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2) );

return beta;
}

我们通过caclNWeights函数计算非终端顶点的权重值，计算公式依据V函数，权重结果放在四个矩阵leftW, upleftW, upW, uprightW中，最后，我们根据像素和权重值构建图，并用max-flow算法解得min-cut，求解的结果放在mask图像中，前景部分的值为GC_PR_FGD，背景部分的值为GC_PR_BGD

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 );
}