opencv最大类间方差法（otsu)实现（三种算法比较）

OTSU算法：就是计算出灰度图最佳阈值的算法

1.先对灰度图进行直方图计算并归一化处理，得到0-255之间每个像素在灰度图中出现的概率，即表示为某个像素在灰度图中出现了n个，灰度图总的像素点为N个，则这个像素的出现概率为Pi=n/N

2.每个灰度图可以由阈值k将灰度图分为A，B两大类，很容易得到A，B类在灰度图中的出现概率以及灰度均值

3.计算灰度图A,B类得类间方差，在最佳阈值K处，求得的类间方差最大，也就是类间方差最大的那个时刻的阈值就为灰度图的最佳阈值

代码实现：

以下三种算法时间复杂度依次降低，前两种算法自己所写，最后一种算法从网上看到大神写的，膜拜！当然如果你想偷懒，可是更直接点用cvThreshold一步实现cvThreshold(srcCopy,thresholdImage,0,255,CV_THRESH_OTSU)，花费时间与算法三相当。

算法一：

处理一张640*480的图片花费时间大概在8000毫秒左右，非常低效，但是是一种非常直接的，显而易见的实现方法。

int otsu(IplImage* src){CvScalar cs = cvSum(src);                             //图像总灰度值int Width = src->width;int Height = src->height;double U_t = 1.0*cs.val[0]/(Width*Height);                  //图像的平均灰度int threshold = 0;double delta = 0;for(int k=0; k<256; k++){double sumOfGray = 0;int num = 0;for(int i=0; i < Width; i++){for(int j=0; j < Height;j++){int t = cvGet2D(src, j,i).val[0];if(t < k){num++;             //统计灰度小于k的像素点的个数sumOfGray += t; //灰度值小于k的像素点的灰度和}}}double w0 = 1.0*num/(Width*Height);                     //灰度值小于阈值k的像素的概率double w1 = 1.0 - w0;                               //灰度值大于k像素的概率double u0 = 1.0*sumOfGray/num;                     //灰度值小于阈值k的像素的平均灰度值double u1 = 1.0*(cs.val[0]-sumOfGray)/(Width*Height-num); //灰度值大于阈值k的像素的平均灰度值double delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //类间方差if(delta_tmp > delta){delta = delta_tmp;threshold = k;}}cout<<cs.val[0]<<endl;cout<<U_t<<endl;return threshold;}

算法二：

借用opencv中的直方图函数，大大减少循环数，同一张图片花费时间大概在16毫秒左右。

int otsu(IplImage* src)  {       int hist_size = 256;    //直方图尺寸        int hist_height = 256;        float range[] = {0,255};  //灰度级的范围        float* ranges[]={range};        //创建一维直方图，统计图像在[0 255]像素的均匀分布        CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);        //计算灰度图像的一维直方图        cvCalcHist(&src,gray_hist,0,0);        //归一化直方图        cvNormalizeHist(gray_hist,1.0);              int Width = src->width;      int Height = src->height;       int threshold = 0;      double delta = 0;CvScalar sumOfGray = cvSum(src);double U_t = sumOfGray.val[0]/(Width*Height);    for(int k=0; k<256; k++)      {          double u0 = 0,  u1 = 0, w0 = 0, w1 = 0, delta_tmp = 0, u00 = 0, u11 = 0;                   for(int i=0; i<k; i++)          {              u00 += cvQueryHistValue_1D(gray_hist,i)*i;   //灰度小于阈值k的像素的平均灰度值              w0 += cvQueryHistValue_1D(gray_hist,i);     //灰度小于阈值k的像素的概率          }u0 = u00/w0;        for(int j=k; j<256; j++)          {              u11 += cvQueryHistValue_1D(gray_hist,j)*j;   //灰度大于阈值k的像素的平均灰度值              w1 += cvQueryHistValue_1D(gray_hist,j);     //灰度小于阈值k的像素的概率          }u1 = u11/w1;        delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2);    //类间方差          if(delta_tmp > delta)          {              delta = delta_tmp;              threshold = k;          }      } cvReleaseHist(&gray_hist);    return threshold;   }

算法三：

一次循环实现，其中的公式需要自己用笔推导一下，非常精妙！完成一次运算的时间毫秒已经不能计量了，显示为0！

int otsu(IplImage* src)  {       int hist_size = 256;    //直方图尺寸        int hist_height = 256;        float range[] = {0,255};  //灰度级的范围        float* ranges[]={range};        //创建一维直方图，统计图像在[0 255]像素的均匀分布        CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);        //计算灰度图像的一维直方图        cvCalcHist(&src,gray_hist,0,0);        //归一化直方图        cvNormalizeHist(gray_hist,1.0);              int Width = src->width;      int Height = src->height;       int threshold = 0;      double delta = 0;double U_t = 0;for(int m=0; m<256; m++){U_t += cvQueryHistValue_1D(gray_hist,m)*m;}double u = 0, w = 0;     for(int k=0; k<256; k++)      {                      u += cvQueryHistValue_1D(gray_hist,k)*k;   //          w += cvQueryHistValue_1D(gray_hist,k);      //灰度大于阈值k的像素的概率  double t = U_t * w - u;double delta_tmp = t * t / (w * (1 - w) );         if(delta_tmp > delta)          {              delta = delta_tmp;              threshold = k;          }      } cvReleaseHist(&gray_hist);    return threshold;   }