opencv:使用高斯混合模型(GMM)源码对视频进行背景差分法

非常感谢thefutureisour对opencv中c++版本的高斯混合模型的源代码完全注释,网上直接使用opencv源码编程的比较少,但是要想自己对高斯混合模型进行优化,或者要想在论文中对高斯混合模型有所创新,必须使用opencv源码来进行编程,而不仅仅是使用opencv的源码接口调用一下修改一下参数。自己废了些脑子提供给网友交流一把

1、 my_background_segm.hpp

#include "opencv2/core/core.hpp"
#include <list>
#include "opencv2/video/background_segm.hpp"    //找到你自己安装包中该文件的位置
namespace cv
{
     class CV_EXPORTS_W my_BackgroundSubtractorMOG : public BackgroundSubtractor
     {
     public:
         //! the default constructor
         CV_WRAP my_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 my_BackgroundSubtractorMOG(int history, int nmixtures, double backgroundRatio, double noiseSigma=0);
         //! the destructor
         virtual ~my_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:
     	//public:
         Size frameSize;
         int frameType;
         Mat bgmodel;
         int nframes;
         int history;
         int nmixtures;
         double varThreshold;
         double backgroundRatio;
         double noiseSigma;
     };
}

2、my_background_segm.cpp

#include "my_background_segm.hpp"
using namespace cv;
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;//默认初始权值
 //不带参数的构造函数,使用默认值
my_BackgroundSubtractorMOG::my_BackgroundSubtractorMOG()
{
    frameSize = Size(0,0);
    frameType = 0;

    nframes = 0;
    nmixtures = defaultNMixtures;
    history = defaultHistory;
    varThreshold = defaultVarThreshold;
    backgroundRatio = defaultBackgroundRatio;
    noiseSigma = defaultNoiseSigma;
}
//带参数的构造函数,使用参数传进来的值
my_BackgroundSubtractorMOG::my_BackgroundSubtractorMOG(int _history, int _nmixtures,
                                                 double _backgroundRatio,
                                                 double _noiseSigma)
{
    frameSize = Size(0,0);
    frameType = 0;

    nframes = 0;
    nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);//不能超过8个,否则就用默认的
    history = _history > 0 ? _history : defaultHistory;//不能小于0,否则就用默认的
    varThreshold = defaultVarThreshold;//门限使用默认的
    backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);//背景门限必须大于0,小于1,否则使用0.95
    noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;//噪声门限大于0
}

my_BackgroundSubtractorMOG::~my_BackgroundSubtractorMOG()
{
}

void my_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 );

    // 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 );//初始化一个1行*XX列的矩阵
	//XX是这样计算的:图像的行*列*混合模型的个数*(1(优先级) + 1(权值) + 2(均值 + 方差) * 通道数)
    bgmodel = Scalar::all(0);//清零
}

//设为模版,就是考虑到了彩色图像与灰度图像两种情况
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 )//如果学习速率为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;//增加权值
						//注意源文章中涉及概率的部分多进行了简化,将概率变为1
                        mptr[k].mean = mu + alpha*diff;//修正均值
                        var = max(var + alpha*(d2 - var), minVar);//开始时方差清零0,所以这里使用噪声方差作为默认方差,否则使用上一次方差
                        mptr[k].var = var;//修正方差
                        mptr[k].sortKey = w/sqrt(var);//重新计算优先级,貌似这里不对,应该使用更新后的mptr[k].weight而不是w

                        for( k1 = k-1; k1 >= 0; k1-- )//从匹配的第k个模型开始向前比较,如果更新后的单高斯模型优先级超过他前面的那个,则交换顺序
                        {
                            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 = min(k, K-1);//有两种情况,当最开始的初始化时,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));//判决:(ucahr)(-true) = 255;(uchar)(-(false)) = 0;
            }
        }
        else//如果学习速率小于等于0,则没有背景更新过程,其他过程类似
        {
            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);
    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(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 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
                {
                    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 my_BackgroundSubtractorMOG::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());//初始化

    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 my_BackgroundSubtractorMOG" );
}

//}

3、main.cpp

#include <iostream>
#include <string>
#include <opencv2/opencv.hpp>
#include "my_background_segm.hpp"   //自己定义的头文件,默认的是直接调用opencv自带的GMM有关的函数,所以本文中重新定义一个不同的类
using namespace cv;
using namespace std;
int count_frame = 0;
int main()
{
  VideoCapture capture("3.avi");
  if( !capture.isOpened() )
  {
    cout<<"读取视频失败"<<endl;
    return -1;
  }
  //获取整个帧数
  long totalFrameNumber = capture.get(CV_CAP_PROP_FRAME_COUNT);
  cout<<"整个视频共"<<totalFrameNumber<<"帧"<<endl;

  //设置开始帧()
  long frameToStart = 1;
  capture.set( CV_CAP_PROP_POS_FRAMES,frameToStart);
  cout<<"从第"<<frameToStart<<"帧开始读"<<endl;

  //设置结束帧
  int frameToStop = 650;

  if(frameToStop < frameToStart)
  {
    cout<<"结束帧小于开始帧,程序错误,即将退出!"<<endl;
    return -1;
  }
  else
  {
    cout<<"结束帧为:第"<<frameToStop<<"帧"<<endl;
  }

  double rate = capture.get(CV_CAP_PROP_FPS);
  int delay = 1000/rate;

  Mat frame;
  //前景图片
  Mat foreground;
  Mat GMM_gray;
  Mat GMM_canny;

  //使用默认参数调用混合高斯模型
  my_BackgroundSubtractorMOG mog;   //使用自己定义的高斯混合模型类
  bool stop(false);
  //currentFrame是在循环体中控制读取到指定的帧后循环结束的变量
  long currentFrame = frameToStart;
  while( !stop )
  {
	  count_frame ++ ;
    if( !capture.read(frame) )
    {
      cout<<"从视频中读取图像失败或者读完整个视频"<<endl;
      return -2;
    }
    cvtColor(frame,GMM_gray,CV_RGB2GRAY);
	//Canny(GMM_gray,GMM_canny,50,150,3);
    //imshow("GMM_canny",GMM_canny);
    imshow("输入视频",frame);
    //参数为:输入图像、输出图像、学习速率
    //mog(GMM_canny,foreground,0.01);
	mog(GMM_gray,foreground,0.01);
	//cout<<mog.nframes<<"  ";
    imshow("前景",foreground);
	medianBlur(foreground,foreground,3);
	imshow("中值滤波后的前景",foreground);

    //按ESC键退出,按其他键会停止在当前帧

    int c = waitKey(delay);

    if ( (char)c == 27 || currentFrame >= frameToStop)
    {
      stop = true;
    }
    if ( c >= 0)
    {
      waitKey(0);
    }
    currentFrame++;

  }

  waitKey(0);
}

据说还用一种是使用cmake对opencv源码进行操作,我没哟使用过,还据说cmake挺简单的,如果有使用的网友可以提出更好的方法哈

时间: 2024-08-30 00:33:47

opencv:使用高斯混合模型(GMM)源码对视频进行背景差分法的相关文章

贝叶斯来理解高斯混合模型GMM

最近学习基础算法<统计学习方法>,看到利用EM算法估计高斯混合模型(GMM)的时候,发现利用贝叶斯的来理解高斯混合模型的应用其实非常合适. 首先,假设我们对于贝叶斯比较熟悉,对高斯分布也熟悉.本文将GMM用于聚类来举例. 除了简单的高斯分布,理论上通过组合多个不同的高斯分布可以构成任意复杂的分布函数.如下图所示: 在最大似然,贝叶斯方法与朴素贝叶斯分类中,2.1中提到高斯概率密度用来计算连续变量情况下的朴素贝叶斯概率.该情况下的高斯分布是训练已知,然后对于输入变量求取其概率密度,结合类别的先验

6. EM算法-高斯混合模型GMM+Lasso详细代码实现

1. 前言 我们之前有介绍过4. EM算法-高斯混合模型GMM详细代码实现,在那片博文里面把GMM说涉及到的过程,可能会遇到的问题,基本讲了.今天我们升级下,主要一起解析下EM算法中GMM(搞事混合模型)带惩罚项的详细代码实现. 2. 原理 由于我们的极大似然公式加上了惩罚项,所以整个推算的过程在几个地方需要修改下. 在带penality的GMM中,我们假设协方差是一个对角矩阵,这样的话,我们计算高斯密度函数的时候,只需要把样本各个维度与对应的\(\mu_k\)和\(\sigma_k\)计算一维

【OpenCV】SIFT原理与源码分析:DoG尺度空间构造

<SIFT原理与源码分析>系列文章索引:http://www.cnblogs.com/tianyalu/p/5467813.html 尺度空间理论 自然界中的物体随着观测尺度不同有不同的表现形态.例如我们形容建筑物用“米”,观测分子.原子等用“纳米”.更形象的例子比如Google地图,滑动鼠标轮可以改变观测地图的尺度,看到的地图绘制也不同:还有电影中的拉伸镜头等等…… 尺度空间中各尺度图像的模糊程度逐渐变大,能够模拟人在距离目标由近到远时目标在视网膜上的形成过程.尺度越大图像越模糊. 为什么要

【OpenCV】SIFT原理与源码分析

SIFT简介 Scale Invariant Feature Transform,尺度不变特征变换匹配算法,是由David G. Lowe在1999年(<Object Recognition from Local Scale-Invariant Features>)提出的高效区域检测算法,在2004年(<Distinctive Image Features from Scale-Invariant Keypoints>)得以完善. SIFT特征对旋转.尺度缩放.亮度变化等保持不变性

OpenCV人脸识别Eigen算法源码分析

1 理论基础 学习Eigen人脸识别算法需要了解一下它用到的几个理论基础,现总结如下: 1.1 协方差矩阵 首先需要了解一下公式: 共公式可以看出:均值描述的是样本集合的平均值,而标准差描述的则是样本集合的各个样本点到均值的距离之平均.以一个国家国民收入为例,均值反映了平均收入,而均方差/方差则反映了贫富差距,如果两个国家国民收入均值相等,则标准差越大说明国家的国民收入越不均衡,贫富差距较大.以上公式都是用来描述一维数据量的,把方差公式推广到二维,则可得到协方差公式: 协方差表明了两个随机变量之

OpenCV人脸识别LBPH算法源码分析

1 背景及理论基础 人脸识别是指将一个需要识别的人脸和人脸库中的某个人脸对应起来(类似于指纹识别),目的是完成识别功能,该术语需要和人脸检测进行区分,人脸检测是在一张图片中把人脸定位出来,完成的是搜寻的功能.从OpenCV2.4开始,加入了新的类FaceRecognizer,该类用于人脸识别,使用它可以方便地进行相关识别实验. 原始的LBP算子定义为在3*3的窗口内,以窗口中心像素为阈值,将相邻的8个像素的灰度值与其进行比较,若周围像素值大于或等于中心像素值,则该像素点的位置被标记为1,否则为0

【OpenCV】SIFT原理与源码分析:关键点描述

<SIFT原理与源码分析>系列文章索引:http://www.cnblogs.com/tianyalu/p/5467813.html 由前一篇<方向赋值>,为找到的关键点即SIFT特征点赋了值,包含位置.尺度和方向的信息.接下来的步骤是关键点描述,即用用一组向量将这个关键点描述出来,这个描述子不但包括关键点,也包括关键点周围对其有贡献的像素点.用来作为目标匹配的依据(所以描述子应该有较高的独特性,以保证匹配率),也可使关键点具有更多的不变特性,如光照变化.3D视点变化等. SIFT

【OpenCV】SIFT原理与源码分析:方向赋值

<SIFT原理与源码分析>系列文章索引:http://www.cnblogs.com/tianyalu/p/5467813.html 由前一篇<关键点搜索与定位>,我们已经找到了关键点.为了实现图像旋转不变性,需要根据检测到的关键点局部图像结构为特征点方向赋值.也就是在findScaleSpaceExtrema()函数里看到的alcOrientationHist()语句: // 计算梯度直方图 float omax = calcOrientationHist(gauss_pyr[o

【OpenCV】SIFT原理与源码分析:关键点搜索与定位

<SIFT原理与源码分析>系列文章索引:http://www.cnblogs.com/tianyalu/p/5467813.html 由前一步<DoG尺度空间构造>,我们得到了DoG高斯差分金字塔: 如上图的金字塔,高斯尺度空间金字塔中每组有五层不同尺度图像,相邻两层相减得到四层DoG结果.关键点搜索就在这四层DoG图像上寻找局部极值点. DoG局部极值点 寻找DoG极值点时,每一个像素点和它所有的相邻点比较,当其大于(或小于)它的图像域和尺度域的所有相邻点时,即为极值点.如下图所