Background Modeling and Foreground Detection -- GMM

Paper

Improved Adaptive Gaussian Mixture Model for Background Subtraction

根据作者Zoran Zivkovic给出的代码修改了一下

Abstract

Background subtraction is a common computer vision task. We analyze the usual pixel-level approach. We develop an efcient adaptive algorithm using Gaussian mixture probability density. Recursive equations are used to constantly update the parameters and but also to simultaneously select the appropriate number of components for each pixel

Gaussian Mixture Model

A.

假设每个像素看做是由M个高斯模型生成，

B.

随着时间推移，迭代更新模型







其中α 是更新率，比如作者源码中设置 α=0.01

om(t)=1 当在GMM中找到“匹配项”(sample is close to a component)，否则om(t)=0

cf : is a measure of the maximum portion of the datathat can belong to foreground objects without inuencing the background model(eg: cf=0.25)

C.

PAMI第九章看了好几遍也弄明白GMM里的update和EM怎么对应起来，

Code

ZoranZivkovicGMM.h

// Author : zengdong_1991
// Date   : 2016-06-13
// HomePage : http://blog.csdn.net/zengdong_1991 // Email  : 782083852@qq.com

/*!
Paper:
【[2004 ICPR] Improved Adaptive Gaussian Mixture Model for Background Subtraction】
*/

#pragma once
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/features2d/features2d.hpp>
#include "math.h"

#define LEARNINGRATE 0.01
#define INITVARIANCE 36.0
#define BGTHRESHOLD 0.75
#define MAXGMMMODELS 5
#define BIASSIGMA_LOW 9
#define BIASSIGMA_HIGHT 18

#define GMMGAUSSIANFLOAT 6    /*#GMM  中float元素个数*/

typedef struct tagGMMGaussian
{
float variance;
float muR;
float muG;
float muB;
float weight;
float significants;    //(weight/sigma)
}GMM, *PGMM;

class ZoranZivkovicGMM
{
public:
//! full constructor
ZoranZivkovicGMM(
float alpha = LEARNINGRATE,
float bgThreshold = BGTHRESHOLD,
float variance = INITVARIANCE,
int maxModels = MAXGMMMODELS,
float biasSigmaLow = BIASSIGMA_LOW,
float biasSigmaHight = BIASSIGMA_HIGHT);
//! default destructor
~ZoranZivkovicGMM();
//!
bool initialize(const cv::Mat& oInitImg, const cv::Mat& oROI);
//!
void operator()(cv::InputArray _image, cv::OutputArray _fgmask, double learningRateOverride = 0);

public:

size_t m_frameNum = 0;
cv::Mat m_oROI;
//! input image size
cv::Size m_oImgSize;
//! input image channel size
size_t m_nImgChannels;
size_t m_height;
size_t m_width;
size_t m_pixelNum;
//! input image type
int m_nImgType;

//列表初始化变量
float m_alpha;
float m_bgThreshold;   //公式（8）
float m_variance;    //sigma0
int m_maxModels;
float m_biasSigmaLow;
float m_biasSigmaHight;

cv::Mat m_modelsPerPixel;        //存储当前各个像素多少个models, 初始化是 0
cv::Mat m_backgroundModels;      //当前背景GMM，  (1, w*h*Num*sizeof(GMM)/sizeof（float）)

};

ZoranZivkovicGMM.cpp

// Author : zengdong_1991
// Date   : 2016-06-13
// HomePage : http://blog.csdn.net/zengdong_1991 // Email  : 782083852@qq.com

#include "ZoranZivkovicGMM.h"
#include <iostream>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

int compareGMM(const void* _gmm1, const void* _gmm2)
{
GMM gmm1 = *(GMM*)_gmm1;
GMM gmm2 = *(GMM*)_gmm2;

if (gmm1.significants < gmm2.significants)
return 1;
else if (gmm1.significants == gmm2.significants)
return 0;
else
return -1;
}

ZoranZivkovicGMM::ZoranZivkovicGMM(float alpha,
float bgThreshold ,
float variance,
int maxModels,
float biasSigmaLow,
float biasSigmaHight)
: m_alpha(alpha)
, m_bgThreshold(bgThreshold)
, m_variance(variance)
, m_maxModels(maxModels)
, m_biasSigmaLow(biasSigmaLow)
, m_biasSigmaHight(biasSigmaHight)
{
std::cout << "#########################Current Parameters Begin###################" << std::endl;
std::cout << "m_alpha = " << m_alpha << std::endl;
std::cout << "m_bgThreshold = " << m_bgThreshold << std::endl;
std::cout << "m_variance = " << m_variance << std::endl;
std::cout << "m_maxModels = " << m_maxModels << std::endl;
std::cout << "m_biasSigmaLow = " << m_biasSigmaLow << std::endl;
std::cout << "m_biasSigmaHight = " << m_biasSigmaHight << std::endl;
std::cout << "#########################Current Parameters End###################" << std::endl;

}

ZoranZivkovicGMM::~ZoranZivkovicGMM()
{

}

bool ZoranZivkovicGMM::initialize(const cv::Mat& oInitImg, const cv::Mat& oROI)
{

CV_Assert(!oInitImg.empty() && oInitImg.cols > 0 && oInitImg.rows > 0);
CV_Assert(oInitImg.isContinuous());
CV_Assert(oInitImg.type() == CV_8UC3);   //必须是彩色图像
m_oImgSize = oInitImg.size();
m_nImgType = oInitImg.type();
m_nImgChannels = oInitImg.channels();
m_height = oInitImg.rows;
m_width = oInitImg.cols;
m_pixelNum = m_height * m_width;

//////////////////////////////////////////////////////////////////////////

//1. m_modelsPerPixel 初始化  表示当前各个像素各有几个GMM
m_modelsPerPixel.create(m_oImgSize, CV_8UC1);
m_modelsPerPixel = cv::Scalar_<uchar>::all(0);

//2. m_backgroundModels 初始化，连续存储GMM模型.
m_backgroundModels.create(1, m_height * m_width * m_maxModels * GMMGAUSSIANFLOAT, CV_32FC1);  //GMMGAUSSIANFLOAT = sizeof(GMM)/sizeof(float)
m_backgroundModels = cv::Scalar_<float>::all(0);

return true;

}

void ZoranZivkovicGMM::operator()(cv::InputArray _image, cv::OutputArray _fgmask, double learningRateOverride)
{

cv::Mat oInputImg = _image.getMat();
CV_Assert(oInputImg.isContinuous());
_fgmask.create(oInputImg.size(), CV_8UC1);
cv::Mat oCurrFGMask = _fgmask.getMat();
oCurrFGMask = cv::Scalar_<uchar>(255);
//用第一帧初始化
bool needToInitialize = m_frameNum == 0 || oInputImg.size() != m_oImgSize || oInputImg.type() != m_nImgType;
if (needToInitialize)
initialize(oInputImg, cv::Mat());
m_frameNum++;

//类型转换
GMM* bgModels = (GMM*)m_backgroundModels.data;

//遍历每个像素
for (size_t x = 0; x < m_height; x++)
{
const uchar* src = oInputImg.ptr<uchar>(x);
uchar* pNumModel = m_modelsPerPixel.ptr<uchar>(x);   //models per pixel
uchar* dstMask = oCurrFGMask.ptr<uchar>(x);

for (size_t y = 0; y < m_width; y++)
{

int posPixel = (x*m_width + y)*m_maxModels;
uchar& numModels = pNumModel[y];   //引用
float sumWeight = 0.0;
cv::Vec3f currPixel(src[y*3], src[y*3 + 1], src[y*3 + 2]);   //当前像素值

//公式（8），B
int backgroundGaussian = 0;
for (size_t i = 0; i < numModels; i++)
{
sumWeight += bgModels[posPixel + i].weight;
if (sumWeight > m_bgThreshold)
{
backgroundGaussian = i + 1;
break;
}
}

bool FitsPDF = false;
sumWeight = 0.0;
for (size_t i = 0; i < numModels; i++)
{
float weight = bgModels[posPixel + i].weight;
if (!FitsPDF)
{

cv::Vec3f bgPixel(bgModels[posPixel + i].muB, bgModels[posPixel + i].muG, bgModels[posPixel + i].muR);
cv::Vec3f diff = currPixel -  bgPixel;
float dist = diff.dot(diff);

if (dist < m_biasSigmaHight * bgModels[posPixel + i].variance && i < backgroundGaussian)
{
dstMask[y] = 0;    //当前像素属于背景
}
if (dist < m_biasSigmaLow * bgModels[posPixel + i].variance)
{
FitsPDF = true;
//按照公式（4）（5）（6）更新参数：
bgModels[posPixel + i].weight += m_alpha*(1 - weight);
bgModels[posPixel + i].muB += (m_alpha / weight)*diff[0];
bgModels[posPixel + i].muG += (m_alpha / weight)*diff[1];
bgModels[posPixel + i].muR += (m_alpha / weight)*diff[2];

float variance = bgModels[posPixel + i].variance;
float sigmaNew = variance + (m_alpha / weight)*(dist - variance);
bgModels[posPixel + i].variance = sigmaNew < 4 ? 4 : sigmaNew > 5 * m_variance ? 5 * m_variance : sigmaNew;
bgModels[posPixel + i].significants = bgModels[posPixel + i].weight / sqrt(bgModels[posPixel + i].variance);
}
else
{
bgModels[posPixel + i].weight *= (1 - m_alpha);
bgModels[posPixel + i].significants = bgModels[posPixel + i].weight / sqrt(bgModels[posPixel + i].variance);
}

}
else
{
bgModels[posPixel + i].weight *= (1 - m_alpha);
bgModels[posPixel + i].significants = bgModels[posPixel + i].weight / sqrt(bgModels[posPixel + i].variance);
}
sumWeight += weight;
}

//weight 归一化
float invSumWeight = 1.0 / sumWeight;
for (size_t i = 0; i < numModels; i++)
{
bgModels[posPixel + i].weight *= invSumWeight;
bgModels[posPixel + i].significants = bgModels[posPixel + i].weight / sqrt(bgModels[posPixel + i].variance);
}

qsort(&bgModels[posPixel], numModels, sizeof(GMM), compareGMM);

if (!FitsPDF)
{
if (numModels < m_maxModels)
{
numModels++;
}

int pos = posPixel + numModels - 1;
bgModels[pos].muB = currPixel[0];
bgModels[pos].muG = currPixel[1];
bgModels[pos].muR = currPixel[2];
bgModels[pos].variance = m_variance;
bgModels[pos].significants = 0;
if (numModels == 1)
{
bgModels[pos].weight = 1.0;
}
else
{
bgModels[pos].weight = m_alpha;
}

//归一化
float sumWeight = 0.0;
for (size_t i = 0; i < numModels; i++)
{
sumWeight += bgModels[posPixel + i].weight;
}
float invSumWeight = 1.0 / sumWeight;
for (size_t i = 0; i < numModels; i++)
{
bgModels[posPixel + i].weight *= invSumWeight;
bgModels[posPixel + i].significants = bgModels[posPixel + i].weight / sqrt(bgModels[posPixel + i].variance);
}
}
qsort(&bgModels[posPixel], numModels, sizeof(GMM), compareGMM);

}

}
}

main.cpp

// Author : zengdong_1991
// Date   : 2016-06-13
// HomePage : http://blog.csdn.net/zengdong_1991 // Email  : 782083852@qq.com

#include <iostream>
#include <cv.h>
#include <highgui.h>

#include "ZoranZivkovicGMM.h"

using namespace std;
using namespace cv;

int main(int argc, char **argv)
{
std::cout << "Using OpenCV " << CV_MAJOR_VERSION << "." << CV_MINOR_VERSION << "." << CV_SUBMINOR_VERSION << std::endl;

CvCapture *capture = 0;
capture = cvCaptureFromAVI("D:\\testVideo\\dataset2014\\highway.avi");  // badminton
if (!capture){
std::cerr << "Cannot open video!" << std::endl;
return 1;
}

ZoranZivkovicGMM bgs;
std::cout << "Press 'q' to quit..." << std::endl;
int key = 0;
IplImage *frame;
bool initialized = false;
cv::Mat img_mask;
cv::Mat img_bkgmodel;

while (key != 'q')
{
frame = cvQueryFrame(capture);
cv::Mat img_input(frame);
bgs(img_input, img_mask);
cv::medianBlur(img_mask, img_mask, 5);
if (!img_mask.empty())
{
cv::medianBlur(img_mask, img_mask, 5);
cv::imshow("output", img_mask);
}
cv::imshow("input", img_input);
cv::waitKey(1);

}
cvDestroyAllWindows();
cvReleaseCapture(&capture);

return 0;
}

BackgroundSubtractorMOG源码(OpenCV)

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                          License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Copyright (C) 2013, OpenCV Foundation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "precomp.hpp"
#include <float.h>

// to make sure we can use these short names
#undef K
#undef L
#undef T

// This is based on the "An Improved Adaptive Background Mixture Model for
// Real-time Tracking with Shadow Detection" by P. KaewTraKulPong and R. Bowden
// http://personal.ee.surrey.ac.uk/Personal/R.Bowden/publications/avbs01/avbs01.pdf //
// The windowing method is used, but not the shadow detection. I make some of my
// own modifications which make more sense. There are some errors in some of their
// equations.
//

namespace cv
{
namespace bgsegm
{

static const int defaultNMixtures = 5;
static const int defaultHistory = 200;
static const double defaultBackgroundRatio = 0.7;
static const double defaultVarThreshold = 2.5*2.5;
static const double defaultNoiseSigma = 30*0.5;
static const double defaultInitialWeight = 0.05;

class BackgroundSubtractorMOGImpl : public BackgroundSubtractorMOG
{
public:
//! the default constructor
BackgroundSubtractorMOGImpl()
{
frameSize = Size(0,0);
frameType = 0;

nframes = 0;
nmixtures = defaultNMixtures;
history = defaultHistory;
varThreshold = defaultVarThreshold;
backgroundRatio = defaultBackgroundRatio;
noiseSigma = defaultNoiseSigma;
name_ = "BackgroundSubtractor.MOG";
}
// the full constructor that takes the length of the history,
// the number of gaussian mixtures, the background ratio parameter and the noise strength
BackgroundSubtractorMOGImpl(int _history, int _nmixtures, double _backgroundRatio, double _noiseSigma=0)
{
frameSize = Size(0,0);
frameType = 0;

nframes = 0;
nmixtures = std::min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);
history = _history > 0 ? _history : defaultHistory;
varThreshold = defaultVarThreshold;
backgroundRatio = std::min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);
noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;
}

//! the update operator
virtual void apply(InputArray image, OutputArray fgmask, double learningRate=0);

//! re-initiaization method
virtual void initialize(Size _frameSize, int _frameType)
{
frameSize = _frameSize;
frameType = _frameType;
nframes = 0;

int nchannels = CV_MAT_CN(frameType);
CV_Assert( CV_MAT_DEPTH(frameType) == CV_8U );

// for each gaussian mixture of each pixel bg model we store ...
// the mixture sort key (w/sum_of_variances), the mixture weight (w),
// the mean (nchannels values) and
// the diagonal covariance matrix (another nchannels values)
bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );
bgmodel = Scalar::all(0);
}

virtual void getBackgroundImage(OutputArray) const
{
CV_Error( Error::StsNotImplemented, "" );
}

virtual int getHistory() const { return history; }
virtual void setHistory(int _nframes) { history = _nframes; }

virtual int getNMixtures() const { return nmixtures; }
virtual void setNMixtures(int nmix) { nmixtures = nmix; }

virtual double getBackgroundRatio() const { return backgroundRatio; }
virtual void setBackgroundRatio(double _backgroundRatio) { backgroundRatio = _backgroundRatio; }

virtual double getNoiseSigma() const { return noiseSigma; }
virtual void setNoiseSigma(double _noiseSigma) { noiseSigma = _noiseSigma; }

virtual void write(FileStorage& fs) const
{
fs << "name" << name_
<< "history" << history
<< "nmixtures" << nmixtures
<< "backgroundRatio" << backgroundRatio
<< "noiseSigma" << noiseSigma;
}

virtual void read(const FileNode& fn)
{
CV_Assert( (String)fn["name"] == name_ );
history = (int)fn["history"];
nmixtures = (int)fn["nmixtures"];
backgroundRatio = (double)fn["backgroundRatio"];
noiseSigma = (double)fn["noiseSigma"];
}

protected:
Size frameSize;
int frameType;
Mat bgmodel;
int nframes;
int history;
int nmixtures;
double varThreshold;
double backgroundRatio;
double noiseSigma;
String name_;
};

template<typename VT> struct MixData
{
float sortKey;
float weight;
VT mean;
VT var;
};

static void process8uC1( const Mat& image, Mat& fgmask, double learningRate,
Mat& bgmodel, int nmixtures, double backgroundRatio,
double varThreshold, double noiseSigma )
{
int x, y, k, k1, rows = image.rows, cols = image.cols;
float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;
int K = nmixtures;
MixData<float>* mptr = (MixData<float>*)bgmodel.data;

const float w0 = (float)defaultInitialWeight;
const float sk0 = (float)(w0/(defaultNoiseSigma*2));
const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);
const float minVar = (float)(noiseSigma*noiseSigma);

for( y = 0; y < rows; y++ )
{
const uchar* src = image.ptr<uchar>(y);
uchar* dst = fgmask.ptr<uchar>(y);

if( alpha > 0 )
{
for( x = 0; x < cols; x++, mptr += K )
{
float wsum = 0;
float pix = src[x];
int kHit = -1, kForeground = -1;

for( k = 0; k < K; k++ )
{
float w = mptr[k].weight;
wsum += w;
if( w < FLT_EPSILON )
break;
float mu = mptr[k].mean;
float var = mptr[k].var;
float diff = pix - mu;
float d2 = diff*diff;
if( d2 < vT*var )
{
wsum -= w;
float dw = alpha*(1.f - w);
mptr[k].weight = w + dw;
mptr[k].mean = mu + alpha*diff;
var = std::max(var + alpha*(d2 - var), minVar);
mptr[k].var = var;
mptr[k].sortKey = w/std::sqrt(var);

for( k1 = k-1; k1 >= 0; k1-- )
{
if( mptr[k1].sortKey >= mptr[k1+1].sortKey )
break;
std::swap( mptr[k1], mptr[k1+1] );
}

kHit = k1+1;
break;
}
}

if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
{
kHit = k = std::min(k, K-1);
wsum += w0 - mptr[k].weight;
mptr[k].weight = w0;
mptr[k].mean = pix;
mptr[k].var = var0;
mptr[k].sortKey = sk0;
}
else
for( ; k < K; k++ )
wsum += mptr[k].weight;

float wscale = 1.f/wsum;
wsum = 0;
for( k = 0; k < K; k++ )
{
wsum += mptr[k].weight *= wscale;
mptr[k].sortKey *= wscale;
if( wsum > T && kForeground < 0 )
kForeground = k+1;
}

dst[x] = (uchar)(-(kHit >= kForeground));
}
}
else
{
for( x = 0; x < cols; x++, mptr += K )
{
float pix = src[x];
int kHit = -1, kForeground = -1;

for( k = 0; k < K; k++ )
{
if( mptr[k].weight < FLT_EPSILON )
break;
float mu = mptr[k].mean;
float var = mptr[k].var;
float diff = pix - mu;
float d2 = diff*diff;
if( d2 < vT*var )
{
kHit = k;
break;
}
}

if( kHit >= 0 )
{
float wsum = 0;
for( k = 0; k < K; k++ )
{
wsum += mptr[k].weight;
if( wsum > T )
{
kForeground = k+1;
break;
}
}
}

dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
}
}
}
}

static void process8uC3( const Mat& image, Mat& fgmask, double learningRate,
Mat& bgmodel, int nmixtures, double backgroundRatio,
double varThreshold, double noiseSigma )
{
int x, y, k, k1, rows = image.rows, cols = image.cols;
float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;
int K = nmixtures;

const float w0 = (float)defaultInitialWeight;
const float sk0 = (float)(w0/(defaultNoiseSigma*2*std::sqrt(3.)));
const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);
const float minVar = (float)(noiseSigma*noiseSigma);
MixData<Vec3f>* mptr = (MixData<Vec3f>*)bgmodel.data;

for( y = 0; y < rows; y++ )
{
const uchar* src = image.ptr<uchar>(y);
uchar* dst = fgmask.ptr<uchar>(y);

if( alpha > 0 )
{
for( x = 0; x < cols; x++, mptr += K )
{
float wsum = 0;
Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
int kHit = -1, kForeground = -1;

for( k = 0; k < K; k++ )
{
float w = mptr[k].weight;
wsum += w;
if( w < FLT_EPSILON )
break;
Vec3f mu = mptr[k].mean;
Vec3f var = mptr[k].var;
Vec3f diff = pix - mu;
float d2 = diff.dot(diff);
if( d2 < vT*(var[0] + var[1] + var[2]) )
{
wsum -= w;
float dw = alpha*(1.f - w);
mptr[k].weight = w + dw;
mptr[k].mean = mu + alpha*diff;
var = Vec3f(std::max(var[0] + alpha*(diff[0]*diff[0] - var[0]), minVar),
std::max(var[1] + alpha*(diff[1]*diff[1] - var[1]), minVar),
std::max(var[2] + alpha*(diff[2]*diff[2] - var[2]), minVar));
mptr[k].var = var;
mptr[k].sortKey = w/std::sqrt(var[0] + var[1] + var[2]);

for( k1 = k-1; k1 >= 0; k1-- )
{
if( mptr[k1].sortKey >= mptr[k1+1].sortKey )
break;
std::swap( mptr[k1], mptr[k1+1] );
}

kHit = k1+1;
break;
}
}

if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
{
kHit = k = std::min(k, K-1);
wsum += w0 - mptr[k].weight;
mptr[k].weight = w0;
mptr[k].mean = pix;
mptr[k].var = Vec3f(var0, var0, var0);
mptr[k].sortKey = sk0;
}
else
for( ; k < K; k++ )
wsum += mptr[k].weight;

float wscale = 1.f/wsum;
wsum = 0;
for( k = 0; k < K; k++ )
{
wsum += mptr[k].weight *= wscale;
mptr[k].sortKey *= wscale;
if( wsum > T && kForeground < 0 )
kForeground = k+1;
}

dst[x] = (uchar)(-(kHit >= kForeground));
}
}
else
{
for( x = 0; x < cols; x++, mptr += K )
{
Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
int kHit = -1, kForeground = -1;

for( k = 0; k < K; k++ )
{
if( mptr[k].weight < FLT_EPSILON )
break;
Vec3f mu = mptr[k].mean;
Vec3f var = mptr[k].var;
Vec3f diff = pix - mu;
float d2 = diff.dot(diff);
if( d2 < vT*(var[0] + var[1] + var[2]) )
{
kHit = k;
break;
}
}

if( kHit >= 0 )
{
float wsum = 0;
for( k = 0; k < K; k++ )
{
wsum += mptr[k].weight;
if( wsum > T )
{
kForeground = k+1;
break;
}
}
}

dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
}
}
}
}

void BackgroundSubtractorMOGImpl::apply(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());

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

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

if( image.type() == CV_8UC1 )
process8uC1( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );
else if( image.type() == CV_8UC3 )
process8uC3( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );
else
CV_Error( Error::StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );
}

Ptr<BackgroundSubtractorMOG> createBackgroundSubtractorMOG(int history, int nmixtures,
double backgroundRatio, double noiseSigma)
{
return makePtr<BackgroundSubtractorMOGImpl>(history, nmixtures, backgroundRatio, noiseSigma);
}

}
}

/* End of file. */