本文详细解析了OpenCV中BackgroundSubtractorMOG类的实现原理与代码细节，介绍了如何利用混合高斯模型进行背景/前景分割，并提供了一种获取背景图像的方法。
本文中的代码选自OpenCV2.4.5的opencv/modules/video/include/opencv2/video/background_segm.hpp和opencv/modules/video/src/bgfg_gaussmix.cpp，对应的类是 cv::BackgroundSubtractorMOG。个人感觉作者的代码写得很不错，而且充分利用了OpenCV的数据类型，内存访问方面也比较高效。下面给出这个cpp文件的代码详解。这个代码的一个明显不足之处就是没有给出获取背景图片的函数，本人在最后给出一种实现方法，算是对原有代码的完善吧。
/*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-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., 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*/
#ifndef __OPENCV_BACKGROUND_SEGM_HPP__
#define __OPENCV_BACKGROUND_SEGM_HPP__
#include "opencv2/core/core.hpp"
#include
namespace cv
{
/*!
The Base Class for Background/Foreground Segmentation
The class is only used to define the common interface for
the whole family of background/foreground segmentation algorithms.
*/
class CV_EXPORTS_W BackgroundSubtractor : public Algorithm
{
public:
//! the virtual destructor
virtual ~BackgroundSubtractor();
//! the update operator that takes the next video frame and returns the current foreground mask as 8-bit binary image.
CV_WRAP_AS(apply) virtual void operator()(InputArray image, OutputArray fgmask,
double learningRate=0);
//! computes a background image
virtual void getBackgroundImage(OutputArray backgroundImage) const;
};
/*!
Gaussian Mixture-based Backbround/Foreground Segmentation Algorithm
The class implements the following algorithm:
"An improved adaptive background mixture model for real-time tracking with shadow detection"
P. KadewTraKuPong and R. Bowden,
Proc. 2nd European Workshp on Advanced Video-Based Surveillance Systems, 2001."
http://personal.ee.surrey.ac.uk/Personal/R.Bowden/publications/avbs01/avbs01.pdf
*/
class CV_EXPORTS_W BackgroundSubtractorMOG : public BackgroundSubtractor
{
public:
//! the default constructor
CV_WRAP BackgroundSubtractorMOG();
//! the full constructor that takes the length of the history, the number of gaussian mixtures, the background ratio parameter and the noise strength
CV_WRAP BackgroundSubtractorMOG(int history, int nmixtures, double backgroundRatio, double noiseSigma=0);
//! the destructor
virtual ~BackgroundSubtractorMOG();
//! the update operator
virtual void operator()(InputArray image, OutputArray fgmask, double learningRate=0);
//! re-initiaization method
virtual void initialize(Size frameSize, int frameType);
virtual AlgorithmInfo* info() const;
protected:
Size frameSize; // 处理视频的尺寸
int frameType; // 每帧图片的类型, 仅支持 CV_8UC1 和 CV_8UC3
Mat bgmodel; // 存储所有像素所有高斯分布模型的矩阵
int nframes; // 已处理的帧数, 用于初始化阶段确定学习率
int history; // 1 / history 表示稳定状态的学习率
int nmixtures; // 每个位置的像素用这么多个高斯分布表示
double varThreshold; // 马氏距离阈值的平方
double backgroundRatio; // 分离前背景的权重累积阈值
double noiseSigma; // 高斯分布的初始方差
};
}
#endif
/*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.
//
//
// Intel License Agreement
//
// Copyright (C) 2000, Intel Corporation, 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 Intel Corporation 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
// 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
{
BackgroundSubtractor::~BackgroundSubtractor() {}
void BackgroundSubtractor::operator()(InputArray, OutputArray, double)
{
}
void BackgroundSubtractor::getBackgroundImage(OutputArray) const
{
}
// 每个位置用这么多个高斯分布建模
static const int defaultNMixtures = 5;
// 稳定状态的学习率等于 1 / defaultHistory
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;
BackgroundSubtractorMOG::BackgroundSubtractorMOG()
{
frameSize = Size(0,0);
frameType = 0;
nframes = 0;
nmixtures = defaultNMixtures;
history = defaultHistory;
varThreshold = defaultVarThreshold;
backgroundRatio = defaultBackgroundRatio;
noiseSigma = defaultNoiseSigma;
}
BackgroundSubtractorMOG::BackgroundSubtractorMOG(int _history, int _nmixtures,
double _backgroundRatio,
double _noiseSigma)
{
frameSize = Size(0,0);
frameType = 0;
nframes = 0;
nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);
history = _history > 0 ? _history : defaultHistory;
varThreshold = defaultVarThreshold;
backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);
noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;
}
BackgroundSubtractorMOG::~BackgroundSubtractorMOG()
{
}
void BackgroundSubtractorMOG::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 );
// 解释一下为什么是 2 + 2 * nchannels
// 首先，每个高斯分布要存储一个权重和一个排序的依据, 即下面的 sort key, 所以需要两个 float 型
// 然后, 每个通道要存储均值和方差, 需要两个 float 型, nchannels 个通道就需要 2 * nchannels 个 float 型
// 综上, 每个高斯分布需要的 float 型数据量是 2 + 2 * nchannels
// bgmodel 是个 cv::Mat 类, 创建矩阵的时候, 行数和列数的取值是没有严格的限制的,
// 只要能够容纳 frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels) 这么多个 float 型数据就行
bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );
bgmodel = Scalar::all(0);
}
// 这个模板类型定义的是一个高斯分布需要存储的数据
// 根据处理图片的不同, VT 的类型不同
// 对于单通道图片, VT = float, 对于单通道图片, VT = cv::Vec3f
template 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;
// mptr 是指向高斯分布的结构体的指针, 因为是单通道数据, 所以 MixData 后面尖括号内的 typename = float
// 后面的代码将按照 MixData 类型访问 bgmodel 中的数据
MixData* mptr = (MixData*)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(y);
uchar* dst = fgmask.ptr(y);
// alpha > 0 不仅检测前景, 而且以 alpha 作为学习率更新背景模型
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;
// 虽然算法指定了每个位置的像素采用 K 个高斯分布表示
// 但是, 注意到, 在初始化的时候, 所有 K 个高斯分布的均值, 方差, 权重和排序依据都初始化为零
// 因此, 权重值等于零可以作为当前下标 k 对应的高斯分布是否已经被使用的依据
// 又因为权重使用浮点数表示, 与零值进行判断不宜采用 w == 0 的形式, 所以有了 w < FLT_EPSILON
// 如果 w < FLT_EPSILON 成立, 则表明, 输入像素值 pix 不能和已有的高斯分布匹配,
// 并且当前位置高斯分布仍然不足 K 个, 因此, 随后需要新增高斯分布
// 这种处理方法避免了每个像素还需要额外的一个整型数据记录当前的有效高斯模型数量, 节省了存储空间
// 但是采用的是浮点数的比较, 而不是整型的比较, 所以从运行时间上说, 是更耗时了
if( w < FLT_EPSILON )
break;
float mu = mptr[k].mean;
float var = mptr[k].var;
float diff = pix - mu;
float d2 = diff*diff;
// 计算马氏距离是否小于阈值
// 其实计算的是马氏距离的平方是否小于阈值的平方, 即 d2 / var < vT
// 为了避免除法运算, 把分母移至不等式的右边
if( d2 < vT*var )
{
// 下面更新被匹配的高斯分布
wsum -= w;
float dw = alpha*(1.f - w);
mptr[k].weight = w + dw;
mptr[k].mean = mu + alpha*diff;
// 注意, 方差的取值不能小于 minVar, 否则会出现意想不到的结果
var = max(var + alpha*(d2 - var), minVar);
mptr[k].var = var;
mptr[k].sortKey = w/sqrt(var);
// 由于被匹配的高斯分布权重增加, 原来按照 sortKey 从大到小排列的高斯分布可能会发生改变
// 接下来把匹配的高斯分布朝着下标变小的方向移动, 直至所有高斯分布都按照 sortKey 从大到小排列
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 是重排序之后被匹配的高斯分布的下标
kHit = k1+1;
break;
}
}
// 没有找到匹配的高斯分布
if( kHit < 0 )
{
// 前面 for 循环中没有找到匹配的高斯分布, 那么就要以当前输入像素值 pix 为均值, 初始方差和初始权重新建一个高斯分布
// 并且将这个被修改的高斯分布的下标赋值给 kHit
// 如果当前位置像素在处理前 K 个高斯分布已经用完,
// 那么要修改的显然是下标为 K - 1 的高斯分布, kHit = K - 1
// 如果当前位置像素的高斯分布仍然不足 K 个,
// 那么要修改的高斯分布是前面由于 w < FLT_EPSILON 而导致 for 循环退出的 k, kHit = k
// 所以有下面一行的三连等式
kHit = k = min(k, K-1);
// 第 k 个高斯分布的权重要变成 w0, 但是在前面的 for 循环中, wsum 加上的是原来的权重
// 所以要减去原来的权重, 再加上新的权重 w0, 于是有下面一行的运算
wsum += w0 - mptr[k].weight;
mptr[k].weight = w0;
mptr[k].mean = pix;
mptr[k].var = var0;
mptr[k].sortKey = sk0;
}
else
// 如果前面的 for 循环中找到了匹配的高斯分布, 匹配的高斯分布之后的那些还没有被累加权重,
// 下面就要把这些权重累加完
for( ; k < K; k++ )
wsum += mptr[k].weight;
// wscale 用来做权重归一化, 确保更新完成之后所有高斯分布权重累加之和等于 1
float wscale = 1.f/wsum;
wsum = 0;
for( k = 0; k < K; k++ )
{
// 权重归一化, 并且重新计算 wsum
wsum += mptr[k].weight *= wscale;
mptr[k].sortKey *= wscale;
// 下面的操作达到的目的是, 下标大于等于 kForeground 的高斯模型对应的是前景
if( wsum > T && kForeground < 0 )
kForeground = k+1;
}
// 如果当前位置属于前景则输出 255, 属于背景则输出 0
// 显然, 属于前景的条件是 kHit >= kForeground
// 如果该条件成立, 运算结果是个整型数据, 用 16 进制数表示就是 0x00000001, 这个数取负数，即补码，得到 0xffffffff
// 强制转换成 uchar 类型, 即 unsigned char 类型, 实际上是截取了最右方的一个字节, 即 0xff, 即十进制的 255
// 如果该条件不成立, 运算结果是个整型数据, 用 16 进制数表示就是 0x00000000, 这个数取负数还是 0x00000000
// 强制转换成 uchar 类型, 即 unsigned char 类型, 实际上是截取了最右方的一个字节, 即 0x00, 即十进制的 0
// 完全可以写成 dst[x] = (kHit >= kForeground) ? 255 : 0;
dst[x] = (uchar)(-(kHit >= kForeground));
}
}
// alpha <= 0, 只检测前景, 不更新背景模型
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*sqrt(3.)));
const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);
const float minVar = (float)(noiseSigma*noiseSigma);
// mptr 是指向高斯分布的结构体的指针, 因为是三通道数据, 所以 MixData 后面尖括号内的 typename = Vec3f
// 后面的代码将按照 MixData 类型访问 bgmodel 中的数据
MixData* mptr = (MixData*)bgmodel.data;
for( y = 0; y < rows; y++ )
{
const uchar* src = image.ptr(y);
uchar* dst = fgmask.ptr(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);
// 下面的的阈值判决采用的实质是近似的马氏距离
// 严格的马氏距离阈值判决是
// diff[0] * diff[0] / var[0] + diff[1] * diff[1] / var[1] + diff[2] + diff[2] / var[2] < vT
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;
// 注意, 每个通道的方差分别运算, 并且确保不小于 minVar
var = Vec3f(max(var[0] + alpha*(diff[0]*diff[0] - var[0]), minVar),
max(var[1] + alpha*(diff[1]*diff[1] - var[1]), minVar),
max(var[2] + alpha*(diff[2]*diff[2] - var[2]), minVar));
mptr[k].var = var;
mptr[k].sortKey = w/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 )
{
kHit = k = 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 BackgroundSubtractorMOG::operator()(InputArray _image, OutputArray _fgmask, double learningRate)
{
Mat image = _image.getMat();
// 该算法采用 operator() 完成初始化和更新的工作
// 如果 nframes, learningRate, frameSize, frameType 取值不合理, 或者和之前的取值不同, 则重新调用 initialize 函数重置参数
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./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( CV_StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );
}
}
/* End of file. */
void BackgroundSubtractorMOG::getBackground(Mat& backImage)
{
if (frameType == CV_8UC1)
{
backImage.create(frameSize, CV_8UC1);
MixData* mptr = (MixData*)bgmodel.data;
unsigned char* ptr = backImage.data;
int length = frameSize.height * frameSize.width;
for (int i = 0; i < length; i++)
{
ptr[0] = saturate_cast(mptr[0].mean);
ptr += 1;
mptr += nmixtures;
}
}
else if (frameType == CV_8UC3)
{
backImage.create(frameSize, CV_8UC3);
MixData* mptr = (MixData*)bgmodel.data;
unsigned char* ptr = backImage.data;
int length = frameSize.height * frameSize.width;
for (int i = 0; i < length; i++)
{
ptr[0] = saturate_cast(mptr[0].mean[0]);
ptr[1] = saturate_cast(mptr[0].mean[1]);
ptr[2] = saturate_cast(mptr[0].mean[2]);
ptr += 3;
mptr += nmixtures;
}
}
}
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
