2-6简单otsu的背景差分法

//2-6简单otsu的背景差分法，这是摄像头实现的效果最好的
/*
最大类间方差法（otsu）的原理：
       阈值将原图象分成前景，背景两个图象。
       前景：用n1,csum, m1来表示在当前阈值下的前景的点数，质量矩，平均灰度
       后景：用n2, sum-csum, m2来表示在当前阈值下的背景的点数，质量矩，平均灰度
       当取最佳阈值时，背景应该与前景差别最大，关键在于如何选择衡量差别的标准
       而在otsu算法中这个衡量差别的标准就是最大类间方差（英文简称otsu，这也就是这个算法名字的来源）
       在本程序中类间方差用sb表示，最大类间方差用fmax
       关于最大类间方差法（otsu）的性能:
       类间方差法对噪音和目标大小十分敏感，它仅对类间方差为单峰的图像产生较好的分割效果。
       当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好，但是类间方差法是用时最少的。
       最大最大类间方差法（otsu）的公式推导：
       记t为前景与背景的分割阈值，前景点数占图像比例为w0， 平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1。
       则图像的总平均灰度为：u=w0*u0+w1*u1。
       前景和背景图象的方差：g=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)=w0*w1*(u0-u1)*(u0-u1),此公式为方差公式，可参照概率论课本
       上面的g的公式也就是下面程序中的sb的表达式
       当方差g最大时，可以认为此时前景和背景差异最大，也就是此时的灰度是最佳阈值.
*/
#include "stdafx.h"
#include 
   
   
    
    
#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;
}