【opencv源码剖析】背景建模mog2

前言

opencv实现的背景建模方法有很多，早期的opencv版本modules/video/src下有acmmm2003(2.3.1版本)、codebook(2.3.1版本)、gmg(2.4.8版本)、gaussmix、gassmix2方法。这两天下了个opencv3.0beta版本，video模块中只有gaussmix、gassmix2被留了下来，其他的方法都被丢到legacy模块去了，legacy模块一般是放即将废弃或老旧的算法的，这也说明了mog、mog2方法依旧老当益壮。(3.0版本还多了个KNN的背景建模算法，暂时没有多余的时间去研究)

网上有很多相关的高斯背景建模原理介绍，就不再多说，这里主要记录下这两天对mog2代码的理解。

不同建模算法对比

运行速度上，mog2有明显的优势，下图引用自OpenCV中背景建模方法mog2——Adaptive GMM算法小结：

	mog	mog2	gmg
电脑1	26904	14386	25533
电脑2	26947	14578	28834

算法结果参考了官方教程，目测效果mog/mog2应该比gmg好点：

Original Frame(Below image shows the 200th frame of a video)：

Result of BackgroundSubtractorMOG：

Result of BackgroundSubtractorMOG2(Gray color region shows shadow region.)：

Result of BackgroundSubtractorGMG：

MOG2算法结构

看了2.3.1、2.4.8、3.0beta三个版本的mog2源码，差别不大。2.4.8代码相对简洁点，3.0beta增加了很多与算法本身无关的接口，因此算法源代码主要参考的是2.4.8版本。下面是算法接口的使用示例：

//Example usage with as cpp class
BackgroundSubtractorMOG2 bg_model;
//For each new image the model is updates using:
bg_model(img, fgmask);

可以看到，完成mog2运算只需两个行代码，一是算法对象构建，而是执行算法。相关源码如下

1 对象构建

有两个构建函数，一个是默认参数，另一个使用你指定的参数，构建主要任务就是参数配置。

BackgroundSubtractorMOG2::BackgroundSubtractorMOG2()
{
    frameSize = Size(0,0);
    frameType = 0;

    nframes = 0;
    history = defaultHistory2;
    varThreshold = defaultVarThreshold2;
    bShadowDetection = 1;

    nmixtures = defaultNMixtures2;
    backgroundRatio = defaultBackgroundRatio2;
    fVarInit = defaultVarInit2;
    fVarMax  = defaultVarMax2;
    fVarMin = defaultVarMin2;

    varThresholdGen = defaultVarThresholdGen2;
    fCT = defaultfCT2;
    nShadowDetection =  defaultnShadowDetection2;
    fTau = defaultfTau;
}

BackgroundSubtractorMOG2::BackgroundSubtractorMOG2(int _history,  float _varThreshold, bool _bShadowDetection)
{
    frameSize = Size(0,0);
    frameType = 0;

    nframes = 0;
    history = _history > 0 ? _history : defaultHistory2;
    varThreshold = (_varThreshold>0)? _varThreshold : defaultVarThreshold2;
    bShadowDetection = _bShadowDetection;

    nmixtures = defaultNMixtures2;
    backgroundRatio = defaultBackgroundRatio2;
    fVarInit = defaultVarInit2;
    fVarMax  = defaultVarMax2;
    fVarMin = defaultVarMin2;

    varThresholdGen = defaultVarThresholdGen2;
    fCT = defaultfCT2;
    nShadowDetection =  defaultnShadowDetection2;
    fTau = defaultfTau;
}

history：如果不手动设置learningRate，history就被用于计算当前的learningRate，此时history越大，learningRate越小，背景更新越慢
varThreshold[Tb]：方差阈值，用于判断当前像素是前景还是背景

nmixtures：高斯模型个数，默认5个

backgroundRatio[TB]：高斯背景模型权重和阈值，nmixtures个模型按权重排序后，只取模型权重累加值大于backgroundRatio的前几个作为背景模型。也就是说如果该值取得非常小，很可能只使用权重最大的高斯模型作为背景(因为仅一个模型权重就大于backgroundRatio了)

fVarInit：新建高斯模型的方差初始值，默认15

fVarMax：背景更新过程中，用于限制高斯模型方差的最大值，默认20

fVarMin：背景更新过程中，用于限制高斯模型方差的最小值，默认4

varThresholdGen[Tg]：方差阈值，用于是否存在匹配的模型，如果不存在则新建一个

fCT：prune?

fTau： Tau is a threshold on how much darker the shadow can be.Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow

2 算法执行

void BackgroundSubtractorMOG2::operator()(InputArray _image, OutputArray _fgmask, double learningRate)
{
    Mat image = _image.getMat();
    bool needToInitialize = nframes == 0 || learningRate >= 1 || image.size() != frameSize || image.type() != frameType;

    if( needToInitialize )
        initialize(image.size(), image.type());

    _fgmask.create( image.size(), CV_8U );
    Mat fgmask = _fgmask.getMat();

    ++nframes;
    learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./min( 2*nframes, history );
    CV_Assert(learningRate >= 0);

    parallel_for_(Range(0, image.rows),
                  MOG2Invoker(image, fgmask,
                              (GMM*)bgmodel.data,
                              (float*)(bgmodel.data + sizeof(GMM)*nmixtures*image.rows*image.cols),
                              bgmodelUsedModes.data, nmixtures, (float)learningRate,
                              (float)varThreshold,
                              backgroundRatio, varThresholdGen,
                              fVarInit, fVarMin, fVarMax, float(-learningRate*fCT), fTau,
                              bShadowDetection, nShadowDetection));
}

在13行可以看到前面说的history与learningRate的计算式，有3种情况：

1) bg_model(img, fgmask)或bg_model(img, fgmask,-1)。函数输入learningRate为-1，learningRate按history值计算，learningRate=1/min(2*nframes,history)。

2) bg_model(img, fgmask,0)。函数输入learningRate为0，背景模型停止更新。

3) bg_model(img, fgmask,n)。n在0~1之间，背景模型更新速度为n，n越大更新越快，算法内部表现为当前帧参与背景更新的权重越大。

mog2算法主要在MOG2Invoker内实现。

3 mog2背景模型数据结构

1)mog2的背景模型数据分配函数：

void BackgroundSubtractorMOG2::initialize(Size _frameSize, int _frameType)
{
    frameSize = _frameSize;
    frameType = _frameType;
    nframes = 0;

    int nchannels = CV_MAT_CN(frameType);
    CV_Assert( nchannels <= CV_CN_MAX );

    // for each gaussian mixture of each pixel bg model we store ...
    // the mixture weight (w),
    // the mean (nchannels values) and
    // the covariance
    bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + nchannels), CV_32F );
    //make the array for keeping track of the used modes per pixel - all zeros at start
    bgmodelUsedModes.create(frameSize,CV_8U);
    bgmodelUsedModes = Scalar::all(0);
}

主要有两个数据块，bgmodel和bgmodelUsedModes。前者即整个算法维护的背景模型数据，后者用于记录对应像素坐标位置所使用的高斯模型个数。

bgmodel的大小为 frameSize.height*frameSize.width*nmixtures*(2 + nchannels)。

从作者的注释中也可以看出，bgmodel里存放的是些什么东西：高斯模型的权重weight，均值mean，方差covariance。这些即是高斯背景建模的核心。

其中权重和方差以下列结构体表示，刚好2个单位的float大小：

struct GMM
{
    float weight;
    float variance;
};

而每个像素通道对应一个均值，刚好有nchannels个单位的float大小。

但要弄清楚这些数据是怎么排列的，还需要看下一个函数，获取mog2背景结果的函数。

2)mog2的背景结果函数：

void BackgroundSubtractorMOG2::getBackgroundImage(OutputArray backgroundImage) const
{
    int nchannels = CV_MAT_CN(frameType);
    CV_Assert( nchannels == 3 );
    Mat meanBackground(frameSize, CV_8UC3, Scalar::all(0));

    int firstGaussianIdx = 0;
    const GMM* gmm = (GMM*)bgmodel.data;
    const Vec3f* mean = reinterpret_cast<const Vec3f*>(gmm + frameSize.width*frameSize.height*nmixtures);
    for(int row=0; row<meanBackground.rows; row++)
    {
        for(int col=0; col<meanBackground.cols; col++)
        {
            int nmodes = bgmodelUsedModes.at<uchar>(row, col);
            Vec3f meanVal;
            float totalWeight = 0.f;
            for(int gaussianIdx = firstGaussianIdx; gaussianIdx < firstGaussianIdx + nmodes; gaussianIdx++)
            {
                GMM gaussian = gmm[gaussianIdx];
                meanVal += gaussian.weight * mean[gaussianIdx];
                totalWeight += gaussian.weight;

                if(totalWeight > backgroundRatio)
                    break;
            }

            meanVal *= (1.f / totalWeight);
            meanBackground.at<Vec3b>(row, col) = Vec3b(meanVal);
            firstGaussianIdx += nmixtures;
        }
    }

    switch(CV_MAT_CN(frameType))
    {
    case 1:
    {
        vector<Mat> channels;
        split(meanBackground, channels);
        channels[0].copyTo(backgroundImage);
        break;
    }

    case 3:
    {
        meanBackground.copyTo(backgroundImage);
        break;
    }

    default:
        CV_Error(CV_StsUnsupportedFormat, "");
    }
}

}

    const Vec3f* mean = reinterpret_cast<const Vec3f*>(gmm + frameSize.width*frameSize.height*nmixtures);

从上面这行可以看出均值mean整块数据是放在gmm(权重、方差)后面的，终于，我们得到了mog2中背景模型bgmodel的数据分布结构：

MOG2的背景更新过程

待续。。

MOG2使用示例

引用自 \安装目录\opencv2.4.8\sources\samples\cpp\bgfg_segm.cpp

//this is a sample for foreground detection functions
int main(int argc, const char** argv)
{
	help();

	CommandLineParser parser(argc, argv, keys);
	bool useCamera = parser.get<bool>("camera");
	string file = parser.get<string>("file_name");
    VideoCapture cap;
    bool update_bg_model = true;

    if( useCamera )
        cap.open(0);
    else
		cap.open(file.c_str());
	parser.printParams();

    if( !cap.isOpened() )
    {
        printf("can not open camera or video file\n");
        return -1;
    }
    
    namedWindow("image", CV_WINDOW_NORMAL);
    namedWindow("foreground mask", CV_WINDOW_NORMAL);
    namedWindow("foreground image", CV_WINDOW_NORMAL);
    namedWindow("mean background image", CV_WINDOW_NORMAL);

    BackgroundSubtractorMOG2 bg_model;
    Mat img, fgmask, fgimg;

    for(;;)
    {
        cap >> img;
        
        if( img.empty() )
            break;
        
        if( fgimg.empty() )
          fgimg.create(img.size(), img.type());

        //update the model
        bg_model(img, fgmask, update_bg_model ? -1 : 0);

        fgimg = Scalar::all(0);
        img.copyTo(fgimg, fgmask);

        Mat bgimg;
        bg_model.getBackgroundImage(bgimg);

        imshow("image", img);
        imshow("foreground mask", fgmask);
        imshow("foreground image", fgimg);
        if(!bgimg.empty())
          imshow("mean background image", bgimg );

        char k = (char)waitKey(30);
        if( k == 27 ) break;
        if( k == ' ' )
        {
            update_bg_model = !update_bg_model;
            if(update_bg_model)
            	printf("Background update is on\n");
            else
            	printf("Background update is off\n");
        }
    }

    return 0;
}

关于blobtracker前景块提取的匀速运动假设

y=a*x+b,这里x是从0~N-1的数值,y分别表示运动轨迹的坐标值x,y,将他们分别进行线性拟合 
根据最小二乘法的线性拟合公式 
a=((x*y)'-x'*y')/((x*x)'-x'*x'); 
b=y'-a*x'; 
这里'是表示平均的意思 

首先有 1^2 + 2^2 + …… + n^2 = n*(n+1)*(2n+1)/6 
有 1 + 2 + …… + n = n*(n+1)/2; 
考虑x，是从0开始则， 
0 + 1^2 + 2^2 + …… + (n-1)^2=n*(n+1)*(2n+1)/6 -n*n=n*(n-1)*(2n-1)/6; 
0 + 1 + 2 + …… + n-1 = n*(n-1)/2; 

公式推导过程 
a = (jsum/N-(N-1)/2*sum/N)/((N*(N+1)*(2*N+1)/6 -N*N)/N-(N-1)/2*(N-1)/2) 
   = (jsum-(N-1)/2*sum)/((N-1)*(2N-1)/6-(N-1)*(N-1)/4)/N 
   = 6*(2*jsum-(N-1)*sum)/((N-1)*(N+1)*N) 
   = 6*(2*jsum-(N-1)*sum)/((N*N-1)*N) 
b = sum/N-a*(N-1)/2; 
   = sum/N-3*(2*jsum-(N-1)*sum)/(N*(N+1)); 
   = ((N+1)*sum-6*jsum-3*(N-1)*sum)/(N*(N+1)); 
   = ((6*N-3)*sum-6*jsum)/(N*(N+1)); 

与程序公式符合 
a[0] = 6*((1-N)*sum[0]+2*jsum[0])/(N*(N*N-1)); 
b[0] = -2*((1-2*N)*sum[0]+3*jsum[0])/(N*(N+1)); 
a[1] = 6*((1-N)*sum[1]+2*jsum[1])/(N*(N*N-1)); 
b[1] = -2*((1-2*N)*sum[1]+3*jsum[1])/(N*(N+1));