OPENCV GrabCut程序详细分析.docx

1、OPENCV GrabCut程序详细分析一、GrabCut函数使用 在OpenCV的源码目录的samples的文件夹下，有grabCut的使用例程，请参考：opencvsamplescppgrabcut.cpp。而grabCut函数的API说明如下：void cv:grabCut( InputArray _img, InputOutputArray _mask, Rect rect, InputOutputArray _bgdModel, InputOutputArray _fgdModel, int iterCount, int mode ) 二、GrabCut源码解读 其中源码包含了gcg

2、raph.hpp这个构建图和max flow/min cut算法的实现文件，这个文件暂时没有解读，后面再更新了。cpp view plaincopyprint? #include precomp.hpp #include gcgraph.hpp #include using namespace cv; class GMM public: static const int componentsCount = 5; GMM( Mat& _model ); double operator()( const Vec3d color ) const; double operator()( int ci,

3、 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 inverseCovscomponentsCount33;

4、/协方差的逆矩阵 double covDetermscomponentsCount; /协方差的行列式 double sumscomponentsCount3; double prodscomponentsCount33; int sampleCountscomponentsCount; int totalSampleCount; ; /背景和前景各有一个对应的GMM（混合高斯模型） GMM:GMM( Mat& _model ) /一个像素的（唯一对应）高斯模型的参数个数或者说一个高斯模型的参数个数 /一个像素RGB三个通道值，故3个均值，3*3个协方差，共用一个权值 const int mo

5、delSize = 3 + 9 + 1; 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_StsBadAr

6、g, _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.ptrdouble(0); /GMM的每个像素的高斯模型的权值变量起始存储指针 mean = coefs + componentsCount; /均值变量起始存储指针 cov = mean + 3*co

7、mponentsCount; /协方差变量起始存储指针 for( int ci = 0; ci 0 ) /计算GMM中第ci个高斯模型的协方差的逆Inverse和行列式Determinant /为了后面计算每个像素属于该高斯模型的概率（也就是数据能量项） calcInverseCovAndDeterm( ci ); /计算一个像素（由color=（B,G,R）三维double型向量来表示）属于这个GMM混合高斯模型的概率。 /也就是把这个像素像素属于componentsCount个高斯模型的概率与对应的权值相乘再相加， /具体见论文的公式（10）。结果从res返回。 /这个相当于计算Gibbs

8、能量的第一个能量项（取负后）。 double GMM:operator()( const Vec3d color ) const double res = 0; for( int ci = 0; ci 0 ) CV_Assert( covDetermsci std:numeric_limitsdouble:epsilon() ); Vec3d diff = color; double* m = mean + 3*ci; diff0 -= m0; diff1 -= m1; diff2 -= m2; double mult = diff0*(diff0*inverseCovsci00 + diff

9、1*inverseCovsci10 + diff2*inverseCovsci20) + diff1*(diff0*inverseCovsci01 + diff1*inverseCovsci11 + diff2*inverseCovsci21) + diff2*(diff0*inverseCovsci02 + diff1*inverseCovsci12 + diff2*inverseCovsci22); res = 1.0f/sqrt(covDetermsci) * exp(-0.5f*mult); return res; /返回这个像素最有可能属于GMM中的哪个高斯模型（概率最大的那个） i

10、nt GMM:whichComponent( const Vec3d color ) const int k = 0; double max = 0; for( int ci = 0; ci max ) k = ci; /找到概率最大的那个，或者说计算结果最大的那个 max = p; return k; /GMM参数学习前的初始化，主要是对要求和的变量置零 void GMM:initLearning() for( int ci = 0; ci componentsCount; ci+) sumsci0 = sumsci1 = sumsci2 = 0; prodsci00 = prodsci01

11、 = prodsci02 = 0; prodsci10 = prodsci11 = prodsci12 = 0; prodsci20 = prodsci21 = prodsci22 = 0; sampleCountsci = 0; totalSampleCount = 0; /增加样本，即为前景或者背景GMM的第ci个高斯模型的像素集（这个像素集是来用估 /计计算这个高斯模型的参数的）增加样本像素。计算加入color这个像素后，像素集 /中所有像素的RGB三个通道的和sums（用来计算均值），还有它的prods（用来计算协方差）， /并且记录这个像素集的像素个数和总的像素个数（用来计算这个高斯

12、模型的权值）。 void GMM:addSample( int ci, const Vec3d color ) sumsci0 += color0; sumsci1 += color1; sumsci2 += color2; prodsci00 += color0*color0; prodsci01 += color0*color1; prodsci02 += color0*color2; prodsci10 += color1*color0; prodsci11 += color1*color1; prodsci12 += color1*color2; prodsci20 += color2

13、*color0; prodsci21 += color2*color1; prodsci22 += color2*color2; sampleCountsci+; totalSampleCount+; /从图像数据中学习GMM的参数：每一个高斯分量的权值、均值和协方差矩阵； /这里相当于论文中“Iterative minimisation”的step 2 void GMM:endLearning() const double variance = 0.01; for( int ci = 0; ci componentsCount; ci+ ) int n = sampleCountsci; /

14、第ci个高斯模型的样本像素个数 if( n = 0 ) coefsci = 0; else /计算第ci个高斯模型的权值系数 coefsci = (double)n/totalSampleCount; /计算第ci个高斯模型的均值 double* m = mean + 3*ci; m0 = sumsci0/n; m1 = sumsci1/n; m2 = sumsci2/n; /计算第ci个高斯模型的协方差 double* c = cov + 9*ci; c0 = prodsci00/n - m0*m0; c1 = prodsci01/n - m0*m1; c2 = prodsci02/n -

15、m0*m2; c3 = prodsci10/n - m1*m0; c4 = prodsci11/n - m1*m1; c5 = prodsci12/n - m1*m2; c6 = prodsci20/n - m2*m0; c7 = prodsci21/n - m2*m1; c8 = prodsci22/n - m2*m2; /计算第ci个高斯模型的协方差的行列式 double dtrm = c0*(c4*c8-c5*c7) - c1*(c3*c8-c5*c6) + c2*(c3*c7-c4*c6); if( dtrm = std:numeric_limitsdouble:epsilon() )

16、 /相当于如果行列式小于等于0，（对角线元素）增加白噪声，避免其变 /为退化（降秩）协方差矩阵（不存在逆矩阵，但后面的计算需要计算逆矩阵）。 / Adds the white noise to avoid singular covariance matrix. c0 += variance; c4 += variance; c8 += variance; /计算第ci个高斯模型的协方差的逆Inverse和行列式Determinant calcInverseCovAndDeterm(ci); /计算协方差的逆Inverse和行列式Determinant void GMM:calcInverseC

17、ovAndDeterm( int ci ) if( coefsci 0 ) /取第ci个高斯模型的协方差的起始指针 double *c = cov + 9*ci; double dtrm = covDetermsci = c0*(c4*c8-c5*c7) - c1*(c3*c8-c5*c6) + c2*(c3*c7-c4*c6); /在C+中，每一种内置的数据类型都拥有不同的属性, 使用库可以获 /得这些基本数据类型的数值属性。因为浮点算法的截断，所以使得，当a=2， /b=3时 10*a/b = 20/b不成立。那怎么办呢？ /这个小正数（epsilon）常量就来了，小正数通常为可用给定数据

18、类型的 /大于1的最小值与1之差来表示。若dtrm结果不大于小正数，那么它几乎为零。 /所以下式保证dtrm0，即行列式的计算正确（协方差对称正定，故行列式大于0）。 CV_Assert( dtrm std:numeric_limitsdouble:epsilon() ); /三阶方阵的求逆 inverseCovsci00 = (c4*c8 - c5*c7) / dtrm; inverseCovsci10 = -(c3*c8 - c5*c6) / dtrm; inverseCovsci20 = (c3*c7 - c4*c6) / dtrm; inverseCovsci01 = -(c1*c8

19、- c2*c7) / dtrm; inverseCovsci11 = (c0*c8 - c2*c6) / dtrm; inverseCovsci21 = -(c0*c7 - c1*c6) / dtrm; inverseCovsci02 = (c1*c5 - c2*c4) / dtrm; inverseCovsci12 = -(c0*c5 - c2*c3) / dtrm; inverseCovsci22 = (c0*c4 - c1*c3) / dtrm; /计算beta，也就是Gibbs能量项中的第二项（平滑项）中的指数项的beta，用来调整 /高或者低对比度时，两个邻域像素的差别的影响的，例如

20、在低对比度时，两个邻域 /像素的差别可能就会比较小，这时候需要乘以一个较大的beta来放大这个差别， /在高对比度时，则需要缩小本身就比较大的差别。 /所以我们需要分析整幅图像的对比度来确定参数beta，具体的见论文公式（5）。 static double calcBeta( const Mat& img ) double beta = 0; for( int y = 0; y img.rows; y+ ) for( int x = 0; x 0 ) / left 0的判断是为了避免在图像边界的时候还计算，导致越界 Vec3d diff = color - (Vec3d)img.at(y,x-

21、1); beta += diff.dot(diff); /矩阵的点乘，也就是各个元素平方的和 if( y0 & x0 ) / upleft Vec3d diff = color - (Vec3d)img.at(y-1,x-1); beta += diff.dot(diff); if( y0 ) / up Vec3d diff = color - (Vec3d)img.at(y-1,x); beta += diff.dot(diff); if( y0 & x/ upright Vec3d diff = color - (Vec3d)img.at(y-1,x+1); beta += diff.do

22、t(diff); if( beta = std:numeric_limitsdouble:epsilon() ) beta = 0; else beta = 1.f / (2 * beta/(4*img.cols*img.rows - 3*img.cols - 3*img.rows + 2) ); /论文公式（5） return beta; /计算图每个非端点顶点（也就是每个像素作为图的一个顶点，不包括源点s和汇点t）与邻域顶点 /的边的权值。由于是无向图，我们计算的是八邻域，那么对于一个顶点，我们计算四个方向就行， /在其他的顶点计算的时候，会把剩余那四个方向的权值计算出来。这样整个图算完后

23、，每个顶点 /与八邻域的顶点的边的权值就都计算出来了。 /这个相当于计算Gibbs能量的第二个能量项（平滑项），具体见论文中公式（4） 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)。 /具

24、体计算时，看下面就明白了 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 i

25、mg.rows; y+ ) for( int x = 0; x =0 ) / left /避免图的边界 Vec3d diff = color - (Vec3d)img.at(y,x-1); leftW.atdouble(y,x) = gamma * exp(-beta*diff.dot(diff); else leftW.atdouble(y,x) = 0; if( x-1=0 & y-1=0 ) / upleft Vec3d diff = color - (Vec3d)img.at(y-1,x-1); upleftW.atdouble(y,x) = gammaDivSqrt2 * exp(-beta*diff.dot(diff); else upleftW.atdouble(y,x) = 0; if( y-1=0 ) / up Vec3d diff = color - (Vec3d)img.at(y-1,x); upW.a