OPENCV GrabCut程序详细分析.docx
《OPENCV GrabCut程序详细分析.docx》由会员分享，可在线阅读，更多相关《OPENCV GrabCut程序详细分析.docx（27页珍藏版）》请在冰豆网上搜索。
OPENCVGrabCut程序详细分析
一、GrabCut函数使用
在OpenCV的源码目录的samples的文件夹下，有grabCut的使用例程，请参考：
opencv\samples\cpp\grabcut.cpp。
而grabCut函数的API说明如下：
voidcv:
:
grabCut（InputArray_img,InputOutputArray_mask,Rectrect,
InputOutputArray_bgdModel,InputOutputArray_fgdModel,
intiterCount,intmode）
二、GrabCut源码解读
其中源码包含了gcgraph.hpp这个构建图和maxflow/mincut算法的实现文件，这个文件暂时没有解读，后面再更新了。
[cpp]viewplaincopyprint?
#include"precomp.hpp"
#include"gcgraph.hpp"
#include
usingnamespacecv;
classGMM
{
public:
staticconstintcomponentsCount=5;
GMM（Mat&_model）;
doubleoperator（）（constVec3dcolor）const;
doubleoperator（）（intci,constVec3dcolor）const;
intwhichComponent（constVec3dcolor）const;
voidinitLearning（）;
voidaddSample（intci,constVec3dcolor）;
voidendLearning（）;
private:
voidcalcInverseCovAndDeterm（intci）;
Matmodel;
double*coefs;
double*mean;
double*cov;
doubleinverseCovs[componentsCount][3][3];//协方差的逆矩阵
doublecovDeterms[componentsCount];//协方差的行列式
doublesums[componentsCount][3];
doubleprods[componentsCount][3][3];
intsampleCounts[componentsCount];
inttotalSampleCount;
};
//背景和前景各有一个对应的GMM（混合高斯模型）
GMM:
:
GMM（Mat&_model）
{
//一个像素的（唯一对应）高斯模型的参数个数或者说一个高斯模型的参数个数
//一个像素RGB三个通道值，故3个均值，3*3个协方差，共用一个权值
constintmodelSize=3+9+1;
if（_model.empty（））
{
//一个GMM共有componentsCount个高斯模型，一个高斯模型有modelSize个模型参数
_model.create（1,modelSize*componentsCount,CV_64FC1）;
_model.setTo（Scalar（0））;
}
elseif（（_model.type（）!
=CV_64FC1）||（_model.rows!
=1）||（_model.cols!
=modelSize*componentsCount））
CV_Error（CV_StsBadArg,"_modelmusthaveCV_64FC1type,rows==1andcols==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（intci=0;ciif（coefs[ci]>0）
//计算GMM中第ci个高斯模型的协方差的逆Inverse和行列式Determinant
//为了后面计算每个像素属于该高斯模型的概率（也就是数据能量项）
calcInverseCovAndDeterm（ci）;
}
//计算一个像素（由color=（B,G,R）三维double型向量来表示）属于这个GMM混合高斯模型的概率。
//也就是把这个像素像素属于componentsCount个高斯模型的概率与对应的权值相乘再相加，
//具体见论文的公式（10）。
结果从res返回。
//这个相当于计算Gibbs能量的第一个能量项（取负后）。
doubleGMM:
:
operator（）（constVec3dcolor）const
{
doubleres=0;
for（intci=0;cires+=coefs[ci]*（*this）（ci,color）;
returnres;
}
//计算一个像素（由color=（B,G,R）三维double型向量来表示）属于第ci个高斯模型的概率。
//具体过程，即高阶的高斯密度模型计算式，具体见论文的公式（10）。
结果从res返回
doubleGMM:
:
operator（）（intci,constVec3dcolor）const
{
doubleres=0;
if（coefs[ci]>0）
{
CV_Assert（covDeterms[ci]>std:
:
numeric_limits<double>:
:
epsilon（））;
Vec3ddiff=color;
double*m=mean+3*ci;
diff[0]-=m[0];diff[1]-=m[1];diff[2]-=m[2];
doublemult=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）;
}
returnres;
}
//返回这个像素最有可能属于GMM中的哪个高斯模型（概率最大的那个）
intGMM:
:
whichComponent（constVec3dcolor）const
{
intk=0;
doublemax=0;
for（intci=0;ci{
doublep=（*this）（ci,color）;
if（p>max）
{
k=ci;//找到概率最大的那个，或者说计算结果最大的那个
max=p;
}
}
returnk;
}
//GMM参数学习前的初始化，主要是对要求和的变量置零
voidGMM:
:
initLearning（）
{
for（intci=0;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（用来计算协方差），
//并且记录这个像素集的像素个数和总的像素个数（用来计算这个高斯模型的权值）。
voidGMM:
:
addSample（intci,constVec3dcolor）
{
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的参数：
每一个高斯分量的权值、均值和协方差矩阵；
//这里相当于论文中“Iterativeminimisation”的step2
voidGMM:
:
endLearning（）
{
constdoublevariance=0.01;
for（intci=0;ci{
intn=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个高斯模型的协方差的行列式
doubledtrm=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，（对角线元素）增加白噪声，避免其变
//为退化（降秩）协方差矩阵（不存在逆矩阵，但后面的计算需要计算逆矩阵）。
//Addsthewhitenoisetoavoidsingularcovariancematrix.
c[0]+=variance;
c[4]+=variance;
c[8]+=variance;
}
//计算第ci个高斯模型的协方差的逆Inverse和行列式Determinant
calcInverseCovAndDeterm（ci）;
}
}
}
//计算协方差的逆Inverse和行列式Determinant
voidGMM:
:
calcInverseCovAndDeterm（intci）
{
if（coefs[ci]>0）
{
//取第ci个高斯模型的协方差的起始指针
double*c=cov+9*ci;
doubledtrm=
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++中，每一种内置的数据类型都拥有不同的属性,使用库可以获
//得这些基本数据类型的数值属性。
因为浮点算法的截断，所以使得，当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）。
staticdoublecalcBeta（constMat&img）
{
doublebeta=0;
for（inty=0;y{
for（intx=0;x{
//计算四个方向邻域两像素的差别，也就是欧式距离或者说二阶范数
//（当所有像素都算完后，就相当于计算八邻域的像素差了）
Vec3dcolor=img.at（y,x）;
if（x>0）//left>0的判断是为了避免在图像边界的时候还计算，导致越界
{
Vec3ddiff=color-（Vec3d）img.at（y,x-1）;
beta+=diff.dot（diff）;//矩阵的点乘，也就是各个元素平方的和
}
if（y>0&&x>0）//upleft
{
Vec3ddiff=color-（Vec3d）img.at（y-1,x-1）;
beta+=diff.dot（diff）;
}
if（y>0）//up
{
Vec3ddiff=color-（Vec3d）img.at（y-1,x）;
beta+=diff.dot（diff）;
}
if（y>0&&x//upright
{
Vec3ddiff=color-（Vec3d）img.at（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））;//论文公式（5）
returnbeta;
}
//计算图每个非端点顶点（也就是每个像素作为图的一个顶点，不包括源点s和汇点t）与邻域顶点
//的边的权值。
由于是无向图，我们计算的是八邻域，那么对于一个顶点，我们计算四个方向就行，
//在其他的顶点计算的时候，会把剩余那四个方向的权值计算出来。
这样整个图算完后，每个顶点
//与八邻域的顶点的边的权值就都计算出来了。
//这个相当于计算Gibbs能量的第二个能量项（平滑项），具体见论文中公式（4）
staticvoidcalcNWeights（constMat&img,Mat&leftW,Mat&upleftW,Mat&upW,
Mat&uprightW,doublebeta,doublegamma）
{
//gammaDivSqrt2相当于公式（4）中的gamma*dis（i,j）^（-1），那么可以知道，
//当i和j是垂直或者水平关系时，dis（i,j）=1，当是对角关系时，dis（i,j）=sqrt（2.0f）。
//具体计算时，看下面就明白了
constdoublegammaDivSqrt2=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（inty=0;y{
for（intx=0;x{
Vec3dcolor=img.at（y,x）;
if（x-1>=0）//left//避免图的边界
{
Vec3ddiff=color-（Vec3d）img.at（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
{
Vec3ddiff=color-（Vec3d）img.at（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
{
Vec3ddiff=color-（Vec3d）img.at（y-1,x）;
upW.a
升级会员 