本文为原创文章，未经本人允许，禁止转载。转载请注明出处。

1.相关资料链接

2.GrabCut算法源代码解读

2.1.载入头文件

#include "precomp.hpp"
#include "opencv2/imgproc/detail/gcgraph.hpp"
#include <limits>

using namespace cv;
using namespace detail;

gcgraph.hpp为最大流最小割算法

2.2.构建GMM类

class GMM
{
public:
    static const int componentsCount = 5;//即高斯混合模型中共有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, double singularFix);
    Mat model;
    
    /*
    每个高斯分量的系数coefs，均值means，协方差cov。
    */
    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;
};

2.3.定义成员函数

2.3.1.定义构造函数

GMM::GMM( Mat& _model )
{
    /*
    一个高斯分量的参数个数，共13个参数：一个像素点有RGB三个通道值，故三个均值，3*3个协方差，每个分量有自己的一个权值。
    */
    const int modelSize = 3/*mean*/ + 9/*covariance*/ + 1/*component weight*/;
    if( _model.empty() )
    {
        _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;
    
    /*
    model中单行参数的排列顺序为：权值、均值、协方差。
    */
    coefs = model.ptr<double>(0);
    mean = coefs + componentsCount;
    cov = mean + 3*componentsCount;

    for( int ci = 0; ci < componentsCount; ci++ )
        if( coefs[ci] > 0 )
             calcInverseCovAndDeterm(ci, 0.0);//计算GMM中第ci个高斯模型的协方差矩阵的逆矩阵和行列式。
    totalSampleCount = 0;
}

opencv中的create函数见下：

Mat::create(int _rows, int _cols, int _type)

用CV_64FC1类型的图存储GMM模型。图只有一行，该行的每个像素值代表一个参数。

如果_model为空，则创建一个_model。行数为1，列数为modelSize*componentsCount，即(3+9+1)*5=65，图像类型为CV_64FC1。并且数值都设为0。

2.3.2.定义operator函数

/*
该像素值属于这个GMM模型（前景或者后景）的概率，即每个高斯分量求得的概率乘相应的权值，然后求和。
*/
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;
}

double GMM::operator()( int ci, const Vec3d color ) const
{
    double res = 0;
    if( coefs[ci] > 0 )//确保每个高斯分量的权值大于0
    {
        /*
        确保各个分量协方差矩阵的行列式的值在一个有效的范围内，否则报错。行列式的值如果等于0，则矩阵不可逆，无法求其逆矩阵。
        */
        CV_Assert( covDeterms[ci] > std::numeric_limits<double>::epsilon() );
        Vec3d diff = color;//diff存放像素点(b,g,r)三个通道的亮度值
        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;
}

假设有:
1
2
double eps = std::numeric_limits<double>::epsilon();
double min = std::numeric_limits<double>::min();
则两者的区别是：
min是编译器认可的最小的正双精度浮点数；
eps是可以保证1.0+eps!=1.0，这个表达式成立的最小的正双精度浮点数。
理论上eps和min两个数字应该一样大不过具体编译器中的定义可能不同。

$double* m$为分量$ci$对应的$B、G、R$的三个均值，即$\mu_b,\mu_g,\mu_r$。

$diff[0]=b-\mu_b;diff[1]=g-\mu_g;diff[2]=r-\mu_r$

$multi=(x-\mu)^T\Sigma^{-1}(x-\mu)$，对应三维高斯模型的指数项的一部分。

$res$计算该像素点属于第$ci$个高斯分量的概率。（省略了$\frac{1}{(2\pi)^{3/2}}$，这样所求得的绝对概率值不正确，但是在求最小化能量函数的时候，常数项可以约掉，不影响最后求的分割结果，还简化了计算）。

2.3.3.确定像素点属于哪个高斯分量

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

某像素点在哪个高斯分量中求得的概率最大，就归为哪个高斯分量。

2.3.4.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;
}

2.3.5.添加样本

void GMM::addSample( int ci, const Vec3d color )
{
    /*
    sum[ci][0]为种子像素点（即硬分割中(前景或者后景)的seeds）的B通道的亮度值的和。
    sum[ci][1]为种子像素点（即硬分割中(前景或者后景)的seeds）的G通道的亮度值的和。
    sum[ci][2]为种子像素点（即硬分割中(前景或者后景)的seeds）的R通道的亮度值的和。
    为求均值做准备。
    */
    sums[ci][0] += color[0]; sums[ci][1] += color[1]; sums[ci][2] += color[2];
    /*
    Cov(X,Y)=E[XY]-E[X]E[Y]。prods求的就相当于ΣXY，为下一步求协方差做准备。
    */
    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];
    /*
    第ci个分量纳入的种子像素点（前景和背景分开计算）的个数。
    用于计算每个高斯分量的权值。
    */
    sampleCounts[ci]++;
    totalSampleCount++;
}

2.3.6.GMM的参数学习

void GMM::endLearning()
{
    CV_Assert(totalSampleCount > 0);//标记的前景seeds和背景seeds的数目都需要大于0。否则报错。
    for( int ci = 0; ci < componentsCount; ci++ )
    {
        /*
        计算每个高斯分量的权值，即属于该分量的像素点数目除以（属于前景或背景的）总的像素点数目。
        */
        int n = sampleCounts[ci];
        if( n == 0 )
            coefs[ci] = 0;
        else
        {
            double inv_n = 1.0 / n;
            coefs[ci] = (double)n/totalSampleCount;

            double* m = mean + 3*ci;
            m[0] = sums[ci][0] * inv_n; m[1] = sums[ci][1] * inv_n; m[2] = sums[ci][2] * inv_n;
            
            /*
            计算第ci个分量对应的协方差矩阵。公式：Cov(X,Y)=E[XY]-E[X]E[Y]。
            */
            double* c = cov + 9*ci;
            c[0] = prods[ci][0][0] * inv_n - m[0]*m[0]; c[1] = prods[ci][0][1] * inv_n - m[0]*m[1]; c[2] = prods[ci][0][2] * inv_n - m[0]*m[2];
            c[3] = prods[ci][1][0] * inv_n - m[1]*m[0]; c[4] = prods[ci][1][1] * inv_n - m[1]*m[1]; c[5] = prods[ci][1][2] * inv_n - m[1]*m[2];
            c[6] = prods[ci][2][0] * inv_n - m[2]*m[0]; c[7] = prods[ci][2][1] * inv_n - m[2]*m[1]; c[8] = prods[ci][2][2] * inv_n - m[2]*m[2];

            calcInverseCovAndDeterm(ci, 0.01);//计算第ci个高斯分量的协方差矩阵的逆矩阵。
        }
    }
}

$m[0],m[1],m[2]为\mu_b,\mu_g,\mu_r$

2.3.7.求协方差矩阵的逆矩阵

void GMM::calcInverseCovAndDeterm(int ci, const double singularFix)
{
    if( coefs[ci] > 0 )
    {
        double *c = cov + 9*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]);//计算第ci个分量的协方差矩阵对应的行列式dtrm。
        if (dtrm <= 1e-6 && singularFix > 0)
        {
            // Adds the white noise to avoid singular covariance matrix.
            c[0] += singularFix;
            c[4] += singularFix;
            c[8] += singularFix;
            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]);
        }
        covDeterms[ci] = dtrm;

        CV_Assert( dtrm > std::numeric_limits<double>::epsilon() );//确保行列式的值在有效范围内
        double inv_dtrm = 1.0 / dtrm;
        inverseCovs[ci][0][0] =  (c[4]*c[8] - c[5]*c[7]) * inv_dtrm;
        inverseCovs[ci][1][0] = -(c[3]*c[8] - c[5]*c[6]) * inv_dtrm;
        inverseCovs[ci][2][0] =  (c[3]*c[7] - c[4]*c[6]) * inv_dtrm;
        inverseCovs[ci][0][1] = -(c[1]*c[8] - c[2]*c[7]) * inv_dtrm;
        inverseCovs[ci][1][1] =  (c[0]*c[8] - c[2]*c[6]) * inv_dtrm;
        inverseCovs[ci][2][1] = -(c[0]*c[7] - c[1]*c[6]) * inv_dtrm;
        inverseCovs[ci][0][2] =  (c[1]*c[5] - c[2]*c[4]) * inv_dtrm;
        inverseCovs[ci][1][2] = -(c[0]*c[5] - c[2]*c[3]) * inv_dtrm;
        inverseCovs[ci][2][2] =  (c[0]*c[4] - c[1]*c[3]) * inv_dtrm;
    }
}

求逆矩阵所用的公式为：$A^{-1}=\frac{1}{\mid A \mid}A^*$

2.2部分和2.3部分用于计算区域项。接下来是边界项的计算。

2.4.计算$\beta$

计算边界项中的$\beta$。

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);//向量的点积，也就是各个元素平方的和。相当于sum((zm-zn)^2)。
            }
            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
        /*
        一个像素点有8个邻居节点，所以一张图的总边数为：(4 * img.cols*img.rows - 3 * img.cols - 3 * img.rows + 2))。
        */
        beta = 1.f / (2 * beta/(4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2) );

    return beta;
}

2.5.计算边界项

计算图每个非端点顶点（也就是每个像素作为图的一个顶点，不包括源点$s$和汇点$t$）与邻域顶点的边的权值。由于是无向图，我们计算的是八邻域，那么对于一个顶点，我们计算四个方向就行，在其他的顶点计算的时候，会把剩余那四个方向的权值计算出来。这样整个图算完后，每个顶点与八邻域的顶点的边的权值就都计算出来了。

这相当于计算$Gibbs$能量的第二个能量项（边界项）。

static void calcNWeights( const Mat& img, Mat& leftW, Mat& upleftW, Mat& upW, Mat& uprightW, double beta, double gamma )
{
    /*
    gammaDivSqrt2相当于公式（4）中的gamma * dist(i,j)^(-1)，那么可以知道，
    当i和j是垂直或者水平关系时，dist(i,j)=1，
    当是对角关系时，dist(i,j)=sqrt(2.0f)。
    */
    const double gammaDivSqrt2 = gamma / std::sqrt(2.0f);
    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));
            }
            else
                leftW.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));
            }
            else
                upleftW.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));
            }
            else
                upW.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));
            }
            else
                uprightW.at<double>(y,x) = 0;
        }
    }
}

分别计算$left,upleft,up,upright$四条边的权值，并以像素的形式存储在对应的图$leftW,upleftW,upW,uprightW$中，如果不存在对应的边，则设权值为0。

2.6.检查mask的正确性

检查$mask$的正确性。$mask$通过用户交互或者程序设定，它是和图像大小一样的单通道图，每个像素点只能取

GC_BGD = 0, // background
GC_FGD = 1, // foreground
GC_PR_BGD = 2, // most probably background
GC_PR_FGD = 3, // most probably foreground

四种枚举值，分别表示该像素点（用户或者程序指定）属于背景、前景、可能为背景或者可能为前景像素。

static void checkMask( const Mat& img, const Mat& mask )
{
    if( mask.empty() )
        CV_Error( CV_StsBadArg, "mask is empty" );//mask不能为空
    if( mask.type() != CV_8UC1 )
        CV_Error( CV_StsBadArg, "mask must have CV_8UC1 type" );//mask的类型必须是8UC1，即8位单通道
    if( mask.cols != img.cols || mask.rows != img.rows )
        CV_Error( CV_StsBadArg, "mask must have as many rows and cols as img" );//mask的大小必须和原图一样
    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 equal "
                    "GC_BGD or GC_FGD or GC_PR_BGD or GC_PR_FGD" );
        }
    }
}

2.7.初始化mask

通过矩形来初始化mask，rect外的全部作为背景，设置为GC_BGD，rect内的设置为GC_PR_FGD（可能为前景）。

static void initMaskWithRect( Mat& mask, Size imgSize, Rect rect )
{
    mask.create( imgSize, CV_8UC1 );
    mask.setTo( GC_BGD );
    
    /*
    建立rect，防止rect的范围越界。
    */
    rect.x = std::max(0, rect.x);
    rect.y = std::max(0, rect.y);
    rect.width = std::min(rect.width, imgSize.width-rect.x);
    rect.height = std::min(rect.height, imgSize.height-rect.y);

    (mask(rect)).setTo( Scalar(GC_PR_FGD) );
}

2.8.初始化GMM

对应论文中“初始化”的步骤

在进行算法第一次迭代时，通过$k-means$算法将GC_BGD区域（或GC_PR_FGD区域）的像素点聚为5类，从而初始化背景GMM和前景GMM模型。

static void initGMMs( const Mat& img, const Mat& mask, GMM& bgdGMM, GMM& fgdGMM )
{
    const int kMeansItCount = 10;
    const int kMeansType = KMEANS_PP_CENTERS;

    Mat bgdLabels, fgdLabels;//建立两个图，分别用于存储前景和背景标签。记录背景和前景的像素样本集中每个像素对应GMM的哪个高斯分量。
    std::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++ )
        {
            /*
            如果该像素点在mask中对应位置上为GC_BGD或者GC_PR_BGD，则将该点纳入到bgdsample中。
            如果该像素点在mask中对应位置上为GC_FGD或者GC_PR_FGD，则将该点纳入到fgdsample中。
            */
            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) );
        }
    }
    CV_Assert( !bgdSamples.empty() && !fgdSamples.empty() );//bgdsample和fgdsample不能为空。
    Mat _bgdSamples( (int)bgdSamples.size(), 3, CV_32FC1, &bgdSamples[0][0] );//构建图_bgdSamples，该图的rows=(int)bgdSamples.size();cols=3;type=CV_32FC1;const Scalar &s=&bgdSamples[0][0]。
    /*
    kmeans++中参数_bgdSamples为：每行一个样本。kmeans的输出为bgdLabels，里面保存的是输入样本集中每一个样本对应的类标签（样本聚为5类后）。
    */
    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 );//同样对属于fgdSamples中的像素点进行聚类。
    
    /*
    基于聚类的结果，构建前景GMM和背景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();
}

这里提供了三种$k-means$算法：

KMEANS_RANDOM_CENTERS=0, // Chooses random centers for k-Means initialization
KMEANS_PP_CENTERS=2, // Uses k-Means++ algorithm for initialization
KMEANS_USE_INITIAL_LABELS=1 // Uses the user-provided labels for K-Means initialization

$OpenCV$的$GrabCut$使用的是k-means++算法。

k-means算法过程如下：

1）从N个文档随机选取K个文档作为质心
2）对剩余的每个文档测量其到每个质心的距离，并把它归到最近的质心的类
3）重新计算已经得到的各个类的质心
4）迭代2～3步直至新的质心与原质心相等或小于指定阈值，算法结束

k-means++算法：
原始k-means算法最开始随机选取数据集中k个点作为聚类中心，而k-means++按照如下的思想选取K个聚类中心：假设已经选取了n个初始聚类中心(0<n<k)，则在选取第n+1个聚类中心时：距离当前n个聚类中心越远的点会有更高的概率被选为第n+1个聚类中心。在选取第一个聚类中心(n=1)时同样通过随机的方法。可以说这也符合我们的直觉：聚类中心当然是互相离得越远越好。这个改进虽然直观简单，但是却非常得有效。

2.9.为每个像素点分配所属的高斯分量

对应论文中“迭代最小化”的第一步

为每个像素分配GMM中所属的高斯模型。

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);
            compIdxs.at<int>(p) = mask.at<uchar>(p) == GC_BGD || mask.at<uchar>(p) == GC_PR_BGD ?
                bgdGMM.whichComponent(color) : fgdGMM.whichComponent(color);
        }
    }
}

如果像素点p=GC_BGD或者p=GC_PR_BGD，则执行$bgdGMM.whichComponent(color)$，返回点$p$最有可能属于背景$GMM$模型中的哪个高斯分量。反之，如果像素点p=GC_FGD或者p=GC_PR_FGD，则执行$fgdGMM.whichComponent(color)$，返回点$p$最有可能属于前景$GMM$模型中的哪个高斯分量。

2.10.GMM参数学习

对应论文中“迭代最小化”的第二步

/*
  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) );
                    else
                        fgdGMM.addSample( ci, img.at<Vec3b>(p) );
                }
            }
        }
    }
    bgdGMM.endLearning();
    fgdGMM.endLearning();
}

2.11.构建GCGraph

通过计算得到的能量项构建图，图的顶点为像素点，图的边由两部分构成，一类边是：每个顶点与汇点t（代表背景）和源点s（代表前景）连接的边，这类边的权值通过Gibbs能量项的第一项能量项来表示。另一类边是：每个顶点与其邻域顶点连接的边，这类边的权值通过Gibbs能量项的第二项能量项来表示。

⚠️构建GCGraph是GrabCut算法中最为耗时的一个步骤。

/*
  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 )
{
    /*
    vtxCount为像素点总数（不包括s点和t点），edgeCount为边的总数（不包括t-link，只有n-link，并且是双向的）。
    */
    int vtxCount = img.cols*img.rows,
        edgeCount = 2*(4*img.cols*img.rows - 3*(img.cols + img.rows) + 2);
    graph.create(vtxCount, edgeCount);//类对象graph（class GCGraph），rows=vtxCount;cols=edgeCount，用于存储点和边的权重。
    Point p;
    for( p.y = 0; p.y < img.rows; p.y++ )
    {
        for( p.x = 0; p.x < img.cols; p.x++)
        {
            // add node
            int vtxIdx = graph.addVtx();//返回点p的索引，按照点加入的顺序，其索引依次为：0,1,2,3,……,vtxCount-1。
            Vec3b color = img.at<Vec3b>(p);

            // set t-weights
            double fromSource, toSink;//用于计算t-link的权重。
            if( mask.at<uchar>(p) == GC_PR_BGD || mask.at<uchar>(p) == GC_PR_FGD )
            {
                /*
                取负对数。这里的log以e为底。
                */
                fromSource = -log( bgdGMM(color) );
                toSink = -log( fgdGMM(color) );
            }
            else if( mask.at<uchar>(p) == GC_BGD )
            {
                fromSource = 0;
                toSink = lambda;
            }
            else // GC_FGD
            {
                fromSource = lambda;
                toSink = 0;
            }//至此，t-link权值计算完毕
            graph.addTermWeights( vtxIdx, fromSource, toSink );

            // set n-weights
            //下面开始设置n-link的权值（之前已经计算好了）。
            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 );
            }
        }
    }
}

2.12.最大流/最小割算法

对应论文中“迭代最小化”的第三步

/*
  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++ )
        {
            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;
                else
                    mask.at<uchar>(p) = GC_PR_BGD;
            }
        }
    }
}

通过图分割的结果来更新$mask$，即最后的图像分割结果。需要注意的是，永远都不会改变用户指定为背景或者前景的像素。

2.13.GrabCut算法API

参数说明：

img——待分割的源图像，必须是8位3通道（CV_8UC3）图像，在处理的过程中不会被修改；
mask——掩码图像，如果使用掩码进行初始化，那么mask保存初始化掩码信息；在执行分割的时候，也可以将用户交互所设定的前景与背景保存到mask中，然后再传入GrabCut函数；在处理结束之后，mask中会保存结果。mask只能取以下四种值：
- GC_BGD（=0），背景；
- GC_FGD（=1），前景；
- GC_PR_BGD（=2），可能的背景；
- GC_PR_FGD（=3），可能的前景。
- 如果没有手工标记GC_BGD或者GC_FGD，那么结果只会有GC_PR_BGD或GC_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 )
{
    CV_INSTRUMENT_REGION();//OpenCV相关算法表现性能测试框架，测量函数执行时间，在函数内部追踪函数执行状况。

    Mat img = _img.getMat();//getMat():把InputArray或者OutputArray转化为Mat
    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 must have CV_8UC3 type" );//输入图像不能为空，类型必须是8UC3。

    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_MASK
            checkMask( img, mask );
        initGMMs( img, mask, bgdGMM, fgdGMM );
    }

    if( iterCount <= 0)
        return;

    if( mode == GC_EVAL_FREEZE_MODEL )
        iterCount = 1;

    if( mode == GC_EVAL || mode == GC_EVAL_FREEZE_MODEL )
        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 );
        if( mode != GC_EVAL_FREEZE_MODEL )
            learnGMMs( img, mask, compIdxs, bgdGMM, fgdGMM );
        constructGCGraph(img, mask, bgdGMM, fgdGMM, lambda, leftW, upleftW, upW, uprightW, graph );
        estimateSegmentation( graph, mask );
    }
}

3.参考资料

1.图像分割之（四）OpenCV的GrabCut函数使用和源码解读

【图像分割】GrabCut算法（源代码解读）

计算机视觉，图像分割，GrabCut，源代码解读