From 0468bdeadde845f718254497dc88d46f900aa4ce Mon Sep 17 00:00:00 2001 From: Vadim Pisarevsky Date: Tue, 28 Dec 2010 16:25:39 +0000 Subject: [PATCH] added background/foreground segmentation algorithm with shadow detection (by Zoran Zivkovic) --- modules/video/src/bgfg_gaussmix2.cpp | 640 +++++++++++++++++++++++++++++++++++ 1 file changed, 640 insertions(+) create mode 100644 modules/video/src/bgfg_gaussmix2.cpp diff --git a/modules/video/src/bgfg_gaussmix2.cpp b/modules/video/src/bgfg_gaussmix2.cpp new file mode 100644 index 0000000..3e6fb85 --- /dev/null +++ b/modules/video/src/bgfg_gaussmix2.cpp @@ -0,0 +1,640 @@ +/*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*/ + +/*//Implementation of the Gaussian mixture model background subtraction from: +// +//"Improved adaptive Gausian mixture model for background subtraction" +//Z.Zivkovic +//International Conference Pattern Recognition, UK, August, 2004 +//http://www.zoranz.net/Publications/zivkovic2004ICPR.pdf +//The code is very fast and performs also shadow detection. +//Number of Gausssian components is adapted per pixel. +// +// and +// +//"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction" +//Z.Zivkovic, F. van der Heijden +//Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006. +// +//The algorithm similar to the standard Stauffer&Grimson algorithm with +//additional selection of the number of the Gaussian components based on: +// +//"Recursive unsupervised learning of finite mixture models " +//Z.Zivkovic, F.van der Heijden +//IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004 +//http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf +// +// +// +//Example usage as part of the CvBGStatModel: +// CvBGStatModel* bg_model = cvCreateGaussianBGModel2( first_frame ); +// +// //update for each frame +// cvUpdateBGStatModel( tmp_frame, bg_model );//segmentation result is in bg_model->foreground +// +// //release at the program termination +// cvReleaseBGStatModel( &bg_model ); +// +//Author: Z.Zivkovic, www.zoranz.net +//Date: 27-April-2005, Version:0.9 +///////////*/ + +#include "cvaux.h" +#include "cvaux_mog2.h" + +int _icvRemoveShadowGMM(long posPixel, + float red, float green, float blue, + unsigned char nModes, + CvPBGMMGaussian* m_aGaussians, + float m_fTb, + float m_fTB, + float m_fTau) +{ + //calculate distances to the modes (+ sort???) + //here we need to go in descending order!!! + long pos; + float tWeight = 0; + float numerator, denominator; + // check all the distributions, marked as background: + for (int iModes=0;iModes= m_fTau))//m_nBeta=1 + { + float dR=a * muR - red; + float dG=a * muG - green; + float dB=a * muB - blue; + + //square distance -slower and less accurate + //float maxDistance = cvSqrt(m_fTb*var); + //if ((fabs(dR) <= maxDistance) && (fabs(dG) <= maxDistance) && (fabs(dB) <= maxDistance)) + //circle + float dist=(dR*dR+dG*dG+dB*dB); + if (dist m_fTB) + { + break; + }; + }; + return 0; +} + +int _icvUpdatePixelBackgroundGMM(long posPixel, + float red, float green, float blue, + unsigned char* pModesUsed, + CvPBGMMGaussian* m_aGaussians, + int m_nM, + float m_fAlphaT, + float m_fTb, + float m_fTB, + float m_fTg, + float m_fSigma, + float m_fPrune) +{ + //calculate distances to the modes (+ sort???) + //here we need to go in descending order!!! + + long pos; + + + bool bFitsPDF=0; + bool bBackground=0; + + float m_fOneMinAlpha=1-m_fAlphaT; + + unsigned char nModes=*pModesUsed; + float totalWeight=0.0f; + + ////// + //go through all modes + for (int iModes=0;iModes20*m_fAlphaT?20*m_fAlphaT:k; + //float sigmanew = var + k*((0.33*(dR*dR+dG*dG+dB*dB))-var); + //float sigmanew = var + k*((dR*dR+dG*dG+dB*dB)-var); + //float sigmanew = var + k*((0.33*dist)-var); + float sigmanew = var + k*(dist-var); + //limit the variance + //m_aGaussians[pos].sigma = sigmanew>70?70:sigmanew; + //m_aGaussians[pos].sigma = sigmanew>5*m_fSigma?5*m_fSigma:sigmanew; + m_aGaussians[pos].sigma =sigmanew< 4 ? 4 : sigmanew>5*m_fSigma?5*m_fSigma:sigmanew; + //m_aGaussians[pos].sigma =sigmanew< 4 ? 4 : sigmanew>3*m_fSigma?3*m_fSigma:sigmanew; + //m_aGaussians[pos].sigma = m_fSigma; + //sort + //all other weights are at the same place and + //only the matched (iModes) is higher -> just find the new place for it + for (int iLocal = iModes;iLocal>0;iLocal--) + { + long posLocal=posPixel + iLocal; + if (weight < (m_aGaussians[posLocal-1].weight)) + { + break; + } + else + { + //swap + CvPBGMMGaussian temp = m_aGaussians[posLocal]; + m_aGaussians[posLocal] = m_aGaussians[posLocal-1]; + m_aGaussians[posLocal-1] = temp; + } + } + + //belongs to the mode + ///// + } + else + { + weight=m_fOneMinAlpha*weight+m_fPrune; + //check prune + if (weight<-m_fPrune) + { + weight=0.0; + nModes--; + // bPrune=1; + //break;//the components are sorted so we can skip the rest + } + } + //check if it fits the current mode (2.5 sigma) + /////// + } + //fit not found yet + ///// + else + { + weight=m_fOneMinAlpha*weight+m_fPrune; + //check prune + if (weight<-m_fPrune) + { + weight=0.0; + nModes--; + } + } + totalWeight+=weight; + m_aGaussians[pos].weight=weight; + } + //go through all modes + ////// + + //renormalize weights + for (int iLocal = 0; iLocal < nModes; iLocal++) + { + m_aGaussians[posPixel+ iLocal].weight = m_aGaussians[posPixel+ iLocal].weight/totalWeight; + } + + //make new mode if needed and exit + if (!bFitsPDF) + { + if (nModes==m_nM) + { + //replace the weakest + } + else + { + //add a new one + nModes++; + } + pos=posPixel+nModes-1; + + if (nModes==1) + m_aGaussians[pos].weight=1; + else + m_aGaussians[pos].weight=m_fAlphaT; + + //renormalize weights + int iLocal; + for (iLocal = 0; iLocal < nModes-1; iLocal++) + { + m_aGaussians[posPixel+ iLocal].weight *=m_fOneMinAlpha; + } + + m_aGaussians[pos].muR=red; + m_aGaussians[pos].muG=green; + m_aGaussians[pos].muB=blue; + m_aGaussians[pos].sigma=m_fSigma; + + //sort + //find the new place for it + for (iLocal = nModes-1;iLocal>0;iLocal--) + { + long posLocal=posPixel + iLocal; + if (m_fAlphaT < (m_aGaussians[posLocal-1].weight)) + { + break; + } + else + { + //swap + CvPBGMMGaussian temp = m_aGaussians[posLocal]; + m_aGaussians[posLocal] = m_aGaussians[posLocal-1]; + m_aGaussians[posLocal-1] = temp; + } + } + } + + //set the number of modes + *pModesUsed=nModes; + + return bBackground; +} + +void _icvReplacePixelBackgroundGMM(long pos, + unsigned char* pData, + CvPBGMMGaussian* m_aGaussians) +{ + pData[0]=(unsigned char) m_aGaussians[pos].muR; + pData[1]=(unsigned char) m_aGaussians[pos].muG; + pData[2]=(unsigned char) m_aGaussians[pos].muB; +} + + +void icvUpdatePixelBackgroundGMM(CvGaussBGStatModel2Data* pGMMData,CvGaussBGStatModel2Params* pGMM, float m_fAlphaT, unsigned char* data,unsigned char* output) +{ + int size=pGMMData->nSize; + unsigned char* pDataCurrent=data; + unsigned char* pUsedModes=pGMMData->rnUsedModes; + unsigned char* pDataOutput=output; + //some constants + int m_nM=pGMM->nM; + //float m_fAlphaT=pGMM->fAlphaT; + + float m_fTb=pGMM->fTb;//Tb - threshold on the Mahalan. dist. + float m_fTB=pGMM->fTB;//1-TF from the paper + float m_fTg=pGMM->fTg;//Tg - when to generate a new component + float m_fSigma=pGMM->fSigma;//initial sigma + float m_fCT=pGMM->fCT;//CT - complexity reduction prior + float m_fPrune=-m_fAlphaT*m_fCT; + float m_fTau=pGMM->fTau; + CvPBGMMGaussian* m_aGaussians=pGMMData->rGMM; + long posPixel=0; + bool m_bShadowDetection=pGMM->bShadowDetection; + unsigned char m_nShadowDetection=pGMM->nShadowDetection; + + //go through the image + for (int i=0;ibRemoveForeground) + { + _icvReplacePixelBackgroundGMM(posPixel,pDataCurrent,m_aGaussians); + } + break; + case 1: + //background + (* pDataOutput)=0; + break; + case 2: + //shadow + (* pDataOutput)=m_nShadowDetection; + if (pGMM->bRemoveForeground) + { + _icvReplacePixelBackgroundGMM(posPixel,pDataCurrent,m_aGaussians); + } + + break; + } + posPixel+=m_nM; + pDataCurrent+=3; + pDataOutput++; + pUsedModes++; + } +} + +////////////////////////////////////////////// +//implementation as part of the CvBGStatModel +static void CV_CDECL icvReleaseGaussianBGModel2( CvGaussBGModel2** bg_model ); +static int CV_CDECL icvUpdateGaussianBGModel2( IplImage* curr_frame, CvGaussBGModel2* bg_model ); + + +CV_IMPL CvBGStatModel* +cvCreateGaussianBGModel2( IplImage* first_frame, CvGaussBGStatModel2Params* parameters ) +{ + CvGaussBGModel2* bg_model = 0; + int w,h,size; + + CV_FUNCNAME( "cvCreateGaussianBGModel2" ); + + __BEGIN__; + + CvGaussBGStatModel2Params params; + + if( !CV_IS_IMAGE(first_frame) ) + CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" ); + + if( !(first_frame->nChannels==3) ) + CV_ERROR( CV_StsBadArg, "Need three channel image (RGB)" ); + + CV_CALL( bg_model = (CvGaussBGModel2*)cvAlloc( sizeof(*bg_model) )); + memset( bg_model, 0, sizeof(*bg_model) ); + bg_model->type = CV_BG_MODEL_MOG2; + bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel2; + bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel2; + + //init parameters + if( parameters == NULL ) + { + /* These constants are defined in cvaux/include/cvaux.h: */ + params.bRemoveForeground=0; + params.bShadowDetection = 1; + params.bPostFiltering=0; + params.minArea=CV_BGFG_MOG2_MINAREA; + + //set parameters + // K - max number of Gaussians per pixel + params.nM = CV_BGFG_MOG2_NGAUSSIANS;//4; + // Tb - the threshold - n var + //pGMM->fTb = 4*4; + params.fTb = CV_BGFG_MOG2_STD_THRESHOLD*CV_BGFG_MOG2_STD_THRESHOLD; + // Tbf - the threshold + //pGMM->fTB = 0.9f;//1-cf from the paper + params.fTB = CV_BGFG_MOG2_BACKGROUND_THRESHOLD; + // Tgenerate - the threshold + params.fTg = CV_BGFG_MOG2_STD_THRESHOLD_GENERATE*CV_BGFG_MOG2_STD_THRESHOLD_GENERATE;//update the mode or generate new + //pGMM->fSigma= 11.0f;//sigma for the new mode + params.fSigma= CV_BGFG_MOG2_SIGMA_INIT; + // alpha - the learning factor + params.fAlphaT=1.0f/CV_BGFG_MOG2_WINDOW_SIZE;//0.003f; + // complexity reduction prior constant + params.fCT=CV_BGFG_MOG2_CT;//0.05f; + + //shadow + // Shadow detection + params.nShadowDetection = CV_BGFG_MOG2_SHADOW_VALUE;//value 0 to turn off + params.fTau = CV_BGFG_MOG2_SHADOW_TAU;//0.5f;// Tau - shadow threshold + } + else + { + params = *parameters; + } + + bg_model->params = params; + + //allocate GMM data + w=first_frame->width; + h=first_frame->height; + size=w*h; + + bg_model->data.nWidth=w; + bg_model->data.nHeight=h; + bg_model->data.nNBands=3; + bg_model->data.nSize=size; + + //GMM for each pixel + bg_model->data.rGMM=(CvPBGMMGaussian*) malloc(size * params.nM * sizeof(CvPBGMMGaussian)); + //used modes per pixel + bg_model->data.rnUsedModes = (unsigned char* ) malloc(size); + memset(bg_model->data.rnUsedModes,0,size);//no modes used + + //prepare storages + CV_CALL( bg_model->background = cvCreateImage(cvSize(w,h), IPL_DEPTH_8U, first_frame->nChannels)); + CV_CALL( bg_model->foreground = cvCreateImage(cvSize(w,h), IPL_DEPTH_8U, 1)); + + //for eventual filtering + CV_CALL( bg_model->storage = cvCreateMemStorage()); + + bg_model->countFrames = 0; + + __END__; + + if( cvGetErrStatus() < 0 ) + { + CvBGStatModel* base_ptr = (CvBGStatModel*)bg_model; + + if( bg_model && bg_model->release ) + bg_model->release( &base_ptr ); + else + cvFree( &bg_model ); + bg_model = 0; + } + + return (CvBGStatModel*)bg_model; +} + + +static void CV_CDECL +icvReleaseGaussianBGModel2( CvGaussBGModel2** _bg_model ) +{ + CV_FUNCNAME( "icvReleaseGaussianBGModel2" ); + + __BEGIN__; + + if( !_bg_model ) + CV_ERROR( CV_StsNullPtr, "" ); + + if( *_bg_model ) + { + CvGaussBGModel2* bg_model = *_bg_model; + + free (bg_model->data.rGMM); + free (bg_model->data.rnUsedModes); + + cvReleaseImage( &bg_model->background ); + cvReleaseImage( &bg_model->foreground ); + cvReleaseMemStorage(&bg_model->storage); + memset( bg_model, 0, sizeof(*bg_model) ); + cvFree( _bg_model ); + } + + __END__; +} + + +static int CV_CDECL +icvUpdateGaussianBGModel2( IplImage* curr_frame, CvGaussBGModel2* bg_model ) +{ + //int i, j, k, n; + int region_count = 0; + CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL; + float alpha,alphaInit; + bg_model->countFrames++; + alpha=bg_model->params.fAlphaT; + + if (bg_model->params.bInit){ + //faster initial updates + alphaInit=(1.0f/(2*bg_model->countFrames+1)); + if (alphaInit>alpha) + { + alpha=alphaInit; + } + else + { + bg_model->params.bInit=0; + } + } + + icvUpdatePixelBackgroundGMM(&bg_model->data,&bg_model->params,alpha,(unsigned char*)curr_frame->imageData,(unsigned char*)bg_model->foreground->imageData); + + if (bg_model->params.bPostFiltering==1) + { + //foreground filtering + + //filter small regions + cvClearMemStorage(bg_model->storage); + + cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 ); + cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 ); + + cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST ); + for( seq = first_seq; seq; seq = seq->h_next ) + { + CvContour* cnt = (CvContour*)seq; + if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea ) + { + //delete small contour + prev_seq = seq->h_prev; + if( prev_seq ) + { + prev_seq->h_next = seq->h_next; + if( seq->h_next ) seq->h_next->h_prev = prev_seq; + } + else + { + first_seq = seq->h_next; + if( seq->h_next ) seq->h_next->h_prev = NULL; + } + } + else + { + region_count++; + } + } + bg_model->foreground_regions = first_seq; + cvZero(bg_model->foreground); + cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1); + + return region_count; + } + else + { + return 1; + } +} + +/* End of file. */ -- 2.7.4