opencv实践程序6——简单（otsu）的背景差分法

//简单otsu的背景差分法，这是摄像头实现的效果最好的
#include <stdio.h>
#include <cv.h>
#include <highgui.h>
using namespace std;
using namespace cv;
int cvOtsu2D(CvMat *pGrayMat)  
{  
    double dHistogram[256][256];        //建立二维灰度直方图  
    double dTrMatrix = 0.0;             //离散矩阵的迹
    int height = pGrayMat->rows;
    int width = pGrayMat->cols;
    int N = height*width;               //总像素数
    int i, j;  
    for(i = 0; i < 256; i++)  
    {  
        for(j = 0; j < 256; j++)  
            dHistogram[i][j] = 0.0;      //初始化变量  
    }  
    for(i = 0; i < height; i++)  
    {  
        for(j = 0; j < width; j++)  
        {  
            unsigned char nData1 = (unsigned char)cvGetReal2D(pGrayMat, i, j);//当前的灰度值  
            unsigned char nData2 = 0;                                                                                                                                
            int nData3 = 0;//注意9个值相加可能超过一个字节                                                                                    
            for(int m = i-1; m <= i+1; m++)  
            {  
                for(int n = j-1; n <= j+1; n++)  
                {  
                    if((m >= 0) && (m < height) && (n >= 0) && (n < width))  
                        nData3 += (unsigned char)cvGetReal2D(pGrayMat, m, n); //当前的灰度值                                                                        
                }  
            }  
            nData2 = (unsigned char)(nData3 / 9);    //对于越界的索引值进行补零,邻域均值  
            dHistogram[nData1][nData2]++;  
        }  
    }  
    for(i = 0; i < 256; i++)  
        for(j = 0; j < 256; j++)  
            dHistogram[i][j] /= N;  //得到归一化的概率分布  
  
    double Pai = 0.0;      //目标区均值矢量i分量  
    double Paj = 0.0;      //目标区均值矢量j分量  
    double Pbi = 0.0;      //背景区均值矢量i分量  
    double Pbj = 0.0;      //背景区均值矢量j分量  
    double Pti = 0.0;      //全局均值矢量i分量  
    double Ptj = 0.0;      //全局均值矢量j分量  
    double W0 = 0.0;       //目标区的联合概率密度  
    double W1 = 0.0;       //背景区的联合概率密度  
    double dData1 = 0.0;  
    double dData2 = 0.0;  
    double dData3 = 0.0;  
    double dData4 = 0.0;   //中间变量  
    int nThreshold_s = 0;  
    int nThreshold_t = 0;  
    double temp = 0.0;     //寻求最大值  
    for(i = 0; i < 256; i++)  
    {  
        for(j = 0; j < 256; j++)  
        {  
            Pti += i*dHistogram[i][j];  
            Ptj += j*dHistogram[i][j];  
        }  
    }  
    for(i = 0; i < 256; i++)  
    {  
        for(j = 0; j < 256; j++)  
        {  
            W0 += dHistogram[i][j];  
            dData1 += i*dHistogram[i][j];  
            dData2 += j*dHistogram[i][j];  
  
            W1 = 1-W0;  
            dData3 = Pti-dData1;  
            dData4 = Ptj-dData2;  
          
            Pai = dData1 / W0;  
            Paj = dData2 / W0;  
            Pbi = dData3 / W1;  
            Pbj = dData4 / W1;   // 得到两个均值向量，用4个分量表示  
            dTrMatrix = ((W0 * Pti - dData1) * (W0 * Pti - dData1) + (W0 * Ptj - dData2) * (W0 * Ptj- dData2)) / (W0 * W1);  
            if(dTrMatrix > temp)  
            {  
                temp = dTrMatrix;  
                nThreshold_s = i;  
                nThreshold_t = j;  
            }  
        }  
    }  
 int   nThreshold = (nThreshold_s + nThreshold_t) / 2;//返回阈值，有多种形式，可以单独某一个，也可                                                         //是两者的均值 
    return nThreshold; 
}  
int main( int argc, char** argv )
{
  IplImage* pFrame = NULL;
  IplImage* pFrImg = NULL;
  IplImage* pBkImg = NULL;


  CvMat* pFrameMat = NULL;
  CvMat* pFrMat = NULL;
  CvMat* pBkMat = NULL;
 
  CvCapture* pCapture = NULL;
 
  int nFrmNum = 0;


  //创建窗口
  cvNamedWindow("video", 1);
  cvNamedWindow("background",1);
  cvNamedWindow("foreground",1);
  //使窗口有序排列
  cvMoveWindow("video", 30, 0);
  cvMoveWindow("background", 360, 0);
  cvMoveWindow("foreground", 690, 0);
   
   pCapture = cvCaptureFromCAM(0);
  // pCapture=cvCaptureFromAVI("2.avi");


    //逐帧读取视频
     while(pFrame = cvQueryFrame( pCapture ))
    {
         nFrmNum++;


         //如果是第一帧，需要申请内存，并初始化
         if(nFrmNum == 1)
        {  
            pBkImg = cvCreateImage(cvSize(pFrame->width, pFrame->height),
                IPL_DEPTH_8U,1);
           pFrImg = cvCreateImage(cvSize(pFrame->width, pFrame->height), 
                IPL_DEPTH_8U,1);
            pBkMat = cvCreateMat(pFrame->height, pFrame->width, CV_32FC1);
           pFrMat = cvCreateMat(pFrame->height, pFrame->width, CV_32FC1);
           pFrameMat = cvCreateMat(pFrame->height, pFrame->width, CV_32FC1);


           //转化成单通道图像再处理
           cvCvtColor(pFrame, pBkImg, CV_BGR2GRAY);
           cvCvtColor(pFrame, pFrImg, CV_BGR2GRAY);


           cvConvert(pFrImg, pFrameMat);
           cvConvert(pFrImg, pFrMat);
           cvConvert(pFrImg, pBkMat);
      }
         else
      {
           cvCvtColor(pFrame, pFrImg, CV_BGR2GRAY);
           cvConvert(pFrImg, pFrameMat);
           //先做高斯滤波，以平滑图像
           //cvSmooth(pFrameMat, pFrameMat, CV_GAUSSIAN, 3, 0, 0);
     
           //当前帧跟背景图相减
           cvAbsDiff(pFrameMat, pBkMat, pFrMat);


           //二值化前景图,cvOtsu2D(pFrMat)
           cvThreshold(pFrMat, pFrImg, cvOtsu2D(pFrMat), 255.0, CV_THRESH_BINARY);
		   


           //更新背景
           cvRunningAvg(pFrameMat, pBkMat, 1, 0);//值设为1对于摄像头来说是效果最好的，但播放视频的时候低一点比较好
           //将背景转化为图像格式，用以显示
           cvConvert(pBkMat, pBkImg);   
                   
           cvShowImage("video", pFrame);
           cvShowImage("background", pBkImg);
           cvShowImage("foreground", pFrImg);


           //如果有按键事件，则跳出循环
           //此等待也为cvShowImage函数提供时间完成显示
           //等待时间可以根据CPU速度调整
         
          if( cvWaitKey(2) >= 0 )
             break;
      }  // end of if-else
    } // end of while-loop


     //销毁窗口
     cvDestroyWindow("video");
     cvDestroyWindow("background");
     cvDestroyWindow("foreground");


     //释放图像和矩阵
     cvReleaseImage(&pFrImg);
     cvReleaseImage(&pBkImg);


     cvReleaseMat(&pFrameMat);
     cvReleaseMat(&pFrMat);
     cvReleaseMat(&pBkMat);


     cvReleaseCapture(&pCapture);
     return 0;
}