OpenCV學習(21) Grabcut算法詳解

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