Background Modeling and Foreground Detection -- SOBS

Paper

1. A Self-Organizing Approach to Background Subtraction for Visual Surveillance Applications
2. The SOBS algorithm: what are the limits?
3. 作者没有公布源码，我按08版Paper简单实现了一下，其中阴影抑制效果太差pass了….

Abstract

We propose an approach based on self organization through artificial neural networks, widely applied in human image processing systems and more generally in cognitive science. The proposed approach can handle scenes containing moving backgrounds, gradual illumination variations and camouflage, has no bootstrapping limitations, can include into the background model shadows cast by moving objects, and achieves robust detection for different types of videos taken with stationary cameras.

Self-Organizing Background Subtraction

The background model constructed and maintained in SOBS algorithm [7] is based on the idea of building the image sequence neural background model by learning in a self-organizing manner image sequence variations, seen as trajectories of pixels in time. The network behaves as a competitive neural network that implements a winner-takeall function, with an associated mechanism that modifies the local synaptic plasticity of neurons, allowing learning to be spatially restricted to the local neighborhood of the most active neurons.

A. Neural Model Representation

给定一输入帧It， 对It中的每一个像素 p 建立一个neural map，这个neural map 是一个n*n的权值向量

比如输入帧大小是N*M，那么最终生成的neural map 的大小是 N*n， M*n
对于输入帧的某个像素p(x,y)，其对应的weight vector 是：

上图的Figure1即是一个(2，3)的图像与其对应的neural map图

B. Neural Model Initialization

Therefore, for each pixel p, the corresponding weight vectors of the model M0(p) are initialized with the pixel brightness value at time t = 0

直接用第一帧初始化，则得到的neural model 可以看成是将第一帧各向row，col扩展n倍得到…

---

C. Background Subtraction and Model Update

At each subsequent time step t, background subtraction is achieved by comparing each pixel p of the t-th sequence frame It with the current pixel model Mt-1(p), in order to determine if there exists a best matching weight vector BM(p) that is close enough to it. If no acceptable matching weight vector exists, p is detected as belonging to a moving object (foreground). Otherwise, if such weight vector is found, it means that p is a background pixel.

1. Finding the Best Match

当前像素It(p)，找到与对应的weight vector中的min（PS：其实不用完全遍历，只要找到了就可以break）

距离度量方式是在HSV空间
——————————————————————————

2. Updating the Model

In order to adapt the background model to scene modifications, including gradual light changes, the model for current pixel p, given by Mt-1(p), should be updated.

如果在weight vector找到满足条件的匹配项 ，假设坐标是（x，y），则对（x，y）周围的值更新

其中

相对于2008年的Paper中，作者在2012年的Paper中提出一种基于Neighborhood Coherence Factor的机制。
具体来说就是2008年的在判断当前像素是前景直接pass，而加入NCF之后就是说当前像素是前景我还要判断邻域的像素属性，如果有50%是背景，那此像素也用于update model….
（PS：Paper里面这个机制实现公式太麻烦了…其实用一个median filter效果是等价的）

——————————————————

3. Parameters

————————————————————————

4. Shadow Detection

阴影抑制
作者使用03年的一篇paper做法

---

---

---

Code

BackgroundSubtractorSOBS.h

// Author : zengdong_1991// Date   : 2016-06-13// HomePage : http://blog.csdn.net/zengdong_1991// Email  : [e-mail]#pragma once#include <opencv2/core/core.hpp>#include <opencv2/imgproc/imgproc.hpp>#include <opencv2/features2d/features2d.hpp>#include "math.h"//! defines the parameters#define WEIGHTS_N   3  //weight vertor 大小 = N^2#define FRAME_NUM_K 45#define EPSILON1 1.0#define EPSILON2 0.008#define PARAM_C1 1.0#define PARAM_C2 0.05#define MAX_W 4.0#define ALPHA_1 (PARAM_C1 / MAX_W)#define ALPHA_2 (PARAM_C2 / MAX_W)#define GAMA 0.7#define BEATA 1.0#define TAU_S 40#define TAU_H 48static const double m_weightW[WEIGHTS_N*WEIGHTS_N] = { 1.0, 2.0, 1.0, 2.0, 4.0, 2.0, 1.0, 2.0, 1.0 };//{  } 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0 class BackgroundSubtractorSOBS {public:    //! full constructor    BackgroundSubtractorSOBS(        int trainFrameNum = FRAME_NUM_K,        double epsilon1 = EPSILON1,        double epsilon2 = EPSILON2,        double c1 = PARAM_C1,        double c2 = PARAM_C2,           double gama = GAMA,        double beta = BEATA,        double tauS = TAU_S,        double tauH = TAU_H,        double maxW = MAX_W);    //! default destructor    virtual ~BackgroundSubtractorSOBS();    //!     virtual void initialize(const cv::Mat& oInitImg, const cv::Mat& oROI);    //!     virtual void operator()(cv::InputArray image, cv::OutputArray fgmask, double learningRateOverride = 0);public:    int m_trainFrameNum;    double m_epsilon1;    double m_epsilon2;    double m_c1;    double m_c2;    double m_maxW;    double m_alpha1;    double m_alpha2;    double m_gama;    double m_beta;    double m_tauS;    double m_tauH;    cv::Mat m_backgroundModel;   //背景模型：  n*rows, n*cols    size_t m_frameNum = 0;    cv::Mat m_ImgHSV;    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;    bool findTheMatchVector(const cv::Vec3b& piexlVal, const size_t x, const size_t y, size_t& offSet_x, size_t& offSet_y, double eta); //true 匹配，false没有匹配    void updateBackground(const cv::Vec3b& piexlVal, const size_t x, const size_t y, size_t& offSet_x, size_t& offSet_y, double alphaT);    double getDistance(const cv::Vec3b& curr, const cv::Vec3b& bg);};

BackgroundSubtractorSOBS.cpp

// Author : zengdong_1991// Date   : 2016-06-13// HomePage : http://blog.csdn.net/zengdong_1991// Email  : [e-mail]#include "BackgroundSubtractorSOBS.h"#include <iostream>#include <opencv2/imgproc/imgproc.hpp>#include <opencv2/highgui/highgui.hpp>BackgroundSubtractorSOBS::BackgroundSubtractorSOBS(int trainFrameNum    ,double epsilon1    ,double epsilon2    ,double c1    ,double c2    ,double gama    ,double beta    ,double tauS    ,double tauH    ,double maxW)    : m_trainFrameNum(trainFrameNum)    , m_epsilon1(epsilon1)    , m_epsilon2(epsilon2)    , m_c1(c1)    , m_c2(c2)    , m_maxW(maxW)    , m_gama(gama)    , m_beta(beta)    , m_tauS(tauS)    , m_tauH(tauH){    m_alpha1 = m_c1 / m_maxW;    m_alpha2 = m_c2 / m_maxW;    std::cout << "#########################Current Parameters Begin###################" << std::endl;    std::cout << "m_trainFrameNum = " << m_trainFrameNum << std::endl;    std::cout << "m_epsilon1 = " << m_epsilon1 << std::endl;    std::cout << "m_epsilon2 = " << m_epsilon2 << std::endl;    std::cout << "m_c1 = " << m_c1 << std::endl;    std::cout << "m_c2 = " << m_c2 << std::endl;    std::cout << "m_maxW = " << m_maxW << std::endl;    std::cout << "m_alpha1 = " << m_alpha1 << std::endl;    std::cout << "m_alpha2 = " << m_alpha2 << std::endl;    std::cout << "m_gama = " << m_gama << std::endl;    std::cout << "m_beta = " << m_beta << std::endl;    std::cout << "m_tauS = " << m_tauS << std::endl;    std::cout << "m_tauH = " << m_tauH << std::endl;    std::cout << "#########################Current Parameters End###################" << std::endl;}BackgroundSubtractorSOBS::~BackgroundSubtractorSOBS() {}void BackgroundSubtractorSOBS::initialize(const cv::Mat& oInitImg, const cv::Mat& oROI){    // == init    CV_Assert(!oInitImg.empty() && oInitImg.cols > 0 && oInitImg.rows > 0);    CV_Assert(oInitImg.isContinuous());    CV_Assert(oInitImg.type() == CV_8UC3);   //必须是彩色图像    if (oInitImg.type() == CV_8UC3)    {        std::vector<cv::Mat> voInitImgChannels;        cv::split(oInitImg, voInitImgChannels);        if (!cv::countNonZero((voInitImgChannels[0] != voInitImgChannels[1]) | (voInitImgChannels[2] != voInitImgChannels[1])))            std::cout << std::endl << "\tBackgroundSubtractorSuBSENSE : Warning, grayscale images should always be passed in CV_8UC1 format for optimal performance." << std::endl;    }    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;    m_frameNum = 1;    //将第一帧转换到 HSV 空间    cv::Mat hsv;    cv::cvtColor(oInitImg, hsv, CV_BGR2HSV);    //H范围是 【0-->180】    //根据hsv初始化背景模型  大小为 (n*rows, n*cols)    m_backgroundModel.create(cv::Size(WEIGHTS_N*m_width, WEIGHTS_N*m_height), CV_8UC3);    for (size_t x = 0; x < m_height; x++)    {        for (size_t y = 0; y < m_width; y++)        {            cv::Vec3b  val = hsv.at<cv::Vec3b>(x, y);            for (size_t i = WEIGHTS_N*x; i < WEIGHTS_N*(x + 1); i++)            {                for (size_t j = WEIGHTS_N*y; j < WEIGHTS_N*(y+1); j++)                {                    m_backgroundModel.at<cv::Vec3b>(i, j) = val;                }            }        }    }    //cv::Mat debugMat = m_backgroundModel;}void BackgroundSubtractorSOBS::operator()(cv::InputArray _image, cv::OutputArray _fgmask, double learningRateOverride){    // == process    cv::Mat oInputImg = _image.getMat();    CV_Assert(oInputImg.type() == m_nImgType && oInputImg.size() == m_oImgSize);    CV_Assert(oInputImg.isContinuous());    _fgmask.create(m_oImgSize, CV_8UC1);    cv::Mat oCurrFGMask = _fgmask.getMat();   //外界 mask 矩阵    oCurrFGMask = cv::Scalar_<uchar>(0);         //转换到HSV空间：     cv::cvtColor(oInputImg, m_ImgHSV, CV_BGR2HSV);    //cv::Mat debugMat = m_ImgHSV;    m_frameNum++;    size_t offSet_x = 0, offSet_y = 0;  //匹配偏移量    double alphaT = m_alpha1 - (m_frameNum*(m_alpha1 - m_alpha2) / m_trainFrameNum); //公式（3）（4）  a(T)    for (size_t x = 0; x < m_height; x++)    {        for (size_t y = 0; y < m_width; y++)        {            cv::Vec3b  val = m_ImgHSV.at<cv::Vec3b>(x, y);            //calibration phase  执行 step 2--> 6            if (m_frameNum <= m_trainFrameNum)            {                //查找匹配                                  if (findTheMatchVector(val, x, y, offSet_x, offSet_y, m_epsilon1))                {                               updateBackground(val, x, y, offSet_x, offSet_y, alphaT);                    oCurrFGMask.at<uchar>(x, y) = 0;                }                else                {                    oCurrFGMask.at<uchar>(x, y) = 255;                }            }            //on line phase  执行 step 2-->10                         else            {                if (findTheMatchVector(val, x, y, offSet_x, offSet_y, m_epsilon2))                {                    updateBackground(val, x, y, offSet_x, offSet_y, m_alpha2);                    oCurrFGMask.at<uchar>(x, y) = 0;                }                /*else if (isShadow(val, x, y, offSet_x, offSet_y))                {                    oCurrFGMask.at<uchar>(x, y) = 0;                        }*/                else                {                    oCurrFGMask.at<uchar>(x, y) = 255;                }            }        }    }}void BackgroundSubtractorSOBS::updateBackground(const cv::Vec3b& piexlVal, const size_t x, const size_t y, size_t& offSet_x, size_t& offSet_y, double alphaT){//  cv::Vec3b* ptrVec = m_backgroundModel.ptr<cv::Vec3b>(0);    int pad = (WEIGHTS_N - 1) / 2;    //当前匹配像素的坐标是    (nx+offset_x, ny+offset_y), 更新其周围的 N*N像素，如下：    for (int i = -pad; i <= pad; i++)    {        int coordinateX = WEIGHTS_N*x + offSet_x + i;        if (coordinateX < 0 || coordinateX >= WEIGHTS_N*m_height)   //边界保护            continue;        for (int j = -pad; j <= pad; j++)        {            int coordinateY = WEIGHTS_N*y + offSet_y + j;            if (coordinateY < 0 || coordinateY >= WEIGHTS_N*m_width)  //边界保护                continue;               cv::Vec3b& pixelValBg = m_backgroundModel.at<cv::Vec3b>(coordinateX, coordinateY);            double alpha = alphaT * m_weightW[(i+pad)*WEIGHTS_N + j + pad];  //【未解决】 w 怎么归一化...现在是【1,2,1；...】            /*防止损失精度还是不这么做好了*///【未解决】 测试以下是否相等...            pixelValBg = (1 - alpha)*pixelValBg + alpha*piexlVal;            //cv::Vec3b At_1 = (1 - alpha)*pixelValBg;            //cv::Vec3b Pt = alpha * piexlVal;            //pixelValBg = At_1 + Pt;            //cv::Vec3b test = pixelValBg;            int a = 0;        }       }}bool BackgroundSubtractorSOBS::findTheMatchVector(const cv::Vec3b& piexlVal, const size_t x, const size_t y, size_t& offSet_x, size_t& offSet_y, double eta){    cv::Mat temp = m_backgroundModel;       //9个像素值...    for (size_t i = 0; i < WEIGHTS_N; i++)    {        for (size_t j = 0; j < WEIGHTS_N; j++)        {            cv::Vec3b piexlValBg = m_backgroundModel.at<cv::Vec3b>(WEIGHTS_N*x + i, WEIGHTS_N*y + j);            //当前像素与背景像素距离小于阈值            double dist = getDistance(piexlVal, piexlValBg);            if (dist < eta)            {                offSet_x = i;                offSet_y = j;                return true;            }        }    }    return false;}double BackgroundSubtractorSOBS::getDistance(const cv::Vec3b& curr, const cv::Vec3b& bg){    //opencv BGR2HSV中的H分量范围属于【0-->180】，将其隐射回【0-->360°(2*pi)】    double hRatio = 2 * CV_PI / 180.0;        double sRatio = 1. / 255.;    double vRatio = 1. / 255.;    double h1 = (double)curr[0] * hRatio;    double s1 = (double)curr[1] * sRatio;    double v1 = (double)curr[2] * vRatio;    double h2 = (double)bg[0] * hRatio;    double s2 = (double)bg[1] * sRatio;    double v2 = (double)bg[2] * vRatio;    double distance = (v1*s1*cos(h1) - v2*s2*cos(h2))*(v1*s1*cos(h1) - v2*s2*cos(h2)) +        (v1*s1*sin(h1) - v2*s2*sin(h2))*(v1*s1*sin(h1) - v2*s2*sin(h2)) +        (v1 - v2)*(v1 - v2);    return (distance);}

main.cpp

// Author : zengdong_1991// Date   : 2016-06-13// HomePage : http://blog.csdn.net/zengdong_1991// Email  : [e-mail]#include <iostream>#include <cv.h>#include <highgui.h>#include "BackgroundSubtractorSOBS.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;    /* Open video file */    CvCapture *capture = 0;    capture = cvCaptureFromAVI("D:\\testVideo\\dataset2014\\highway.avi");    if (!capture){        std::cerr << "Cannot open video!" << std::endl;        return 1;    }    BackgroundSubtractorSOBS 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);        if (!initialized)        {            bgs.initialize(img_input, cv::Mat());            initialized = true;            continue;        }        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);        //    key = cvWaitKey(1);    }    cvDestroyAllWindows();    cvReleaseCapture(&capture);    return 0;}

---

---

---

---

后记

1. 这算法似曾相识啊
2. 教训： 代码优化一定要等到实现完功能之后再做…..
   eg: .at(i, j)