otsu阈值处理

otsu法（最大类间方差法，有时也称之为大津算法）使用的是聚类的思想，把图像的灰度数按灰度级分成2个部分，使得两个部分之间的灰度值差异最大，每个部分之间的灰度差异最小，通过方差的计算来寻找一个合适的灰度级别 来划分。  所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响。因此,使类间方差最大的分割意味着错分概率最小。

设t为设定的阈值。

wo： 分开后  前景像素点数占图像的比例

uo：  分开后  前景像素点的平均灰度

w1：分开后  背景像素点数占图像的比例

u1：  分开后  背景像素点的平均灰度

u=w0*u0 + w1*u1 :图像总平均灰度

从L个灰度级遍历t，使得t为某个值的时候，前景和背景的方差最大， 则 这个 t  值便是我们要求得的阈值。

其中，方差的计算公式如下：

g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)

[             此公式计算量较大，可以采用：      g = wo * w1 * (uo - u1) * (uo - u1)                ]

由于otsu算法是对图像的灰度级进行聚类，所以  在执行otsu算法之前，需要计算该图像的灰度直方图。

#include <stdio.h>  

#include <cv.h>  
#include <highgui.h>  
#include <math.h>  
  
#pragma comment(lib, "cv.lib")  
#pragma comment(lib, "cxcore.lib")  
#pragma comment(lib, "highgui.lib")  
int Otsu(IplImage* src);  
  
int main()  
{  
    IplImage* img = cvLoadImage("00.jpg",0);  
    IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);  
    int threshold = Otsu(img);  
  
    cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);  
  
  
    cvNamedWindow( "img", 1 );  
    cvShowImage("img", dst);  
  
  
    cvWaitKey(-1);  
  
    cvReleaseImage(&img);  
    cvReleaseImage(&dst);  
  
    cvDestroyWindow( "dst" );  
    return 0;  
}  
  
int Otsu(IplImage* src)    
{    
int m,n;
int threshold;      
    long sum0 = 0, sum1 = 0; //存储前景的灰度总和和背景灰度总和  
    long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数  
    double w0 = 0, w1 = 0; //前景和背景所占整幅图像的比例  
    double u0 = 0, u1 = 0;  //前景和背景的平均灰度  
    double variance = 0; //最大类间方差  
    int i, j;  
    double u = 0;  
    double maxVariance = 0; 
    int height=src->height;    
    int width=src->width;        
    long size = height * width;   
  
    //histogram    
    float histogram[256] = {0};    
    for(m=0; m < height; m++)  
    {    
        unsigned char* p=(unsigned char*)src->imageData + src->widthStep * m;    
        for(n = 0; n < width; n++)   
        {    
            histogram[(*p++)]++;    //注：( ++的优先级大于* i=*p++ 相当于 i=*(p++)，又相当于 i=*p; p++;)

        }    
    }    
  
    
    for(i = 1; i < 256; i++) //一次遍历每个像素  
    {    
        sum0 = 0;  
        sum1 = 0;   
        cnt0 = 0;  
        cnt1 = 0;  
        w0 = 0;  
        w1 = 0;  
        for(j = 0; j < i; j++)  
        {  
            cnt0 += histogram[j];  
            sum0 += j * histogram[j];  
        }  
  
        u0 = (double)sum0 /  cnt0;   
        w0 = (double)cnt0 / size;  
  
        for(j = i ; j <= 255; j++)  
        {  
            cnt1 += histogram[j];  
            sum1 += j * histogram[j];  
        }  
  
        u1 = (double)sum1 / cnt1;  
        w1 = 1 - w0; // (double)cnt1 / size;  
  
        u = u0 * w0 + u1 * w1; //图像的平均灰度  
        printf("u = %f\n", u);  
        //variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  
        variance =  w0 * w1 *  (u0 - u1) * (u0 - u1);  
        if(variance > maxVariance)   
        {    
            maxVariance = variance;    
            threshold = i;    
        }   
    }    
  
    printf("threshold = %d\n", threshold);  
    return threshold;    
}