一、Otsu算法原理

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

设t为设定的阈值。

w0	分开后前景像素点数占图像的比例
u0	分开后前景像素点的平均灰度
w1	分开后背景像素点数占图像的比例
u1	分开后背景像素点的平均灰度

 

 

 

 

图像总平均灰度为： u = w0∗u0 + w1∗u1 

 

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

 g = wo∗(u0−u)∗(u0−u) + w1∗(u1−u)∗(u1−u) 

此公式计算量较大，可以采用：

 g = w0∗w1∗(u0−u1)∗(u0−u1) 

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

二、代码实现算法

第一份代码：来自Augusdi

1 #include
     
       2 #include
      
        3 //#include
       
         4 #include
        
          5  6 int Otsu(IplImage* src); 7  8 int main() 9 {10     IplImage* img = cvLoadImage("E:\\house.jpg", 0);11     IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);12     int threshold = Otsu(img);13 14     cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);15 16 17     cvNamedWindow("img", 1);18     cvShowImage("img", dst);19 20 21     cvWaitKey(-1);22 23     cvReleaseImage(&img);24     cvReleaseImage(&dst);25 26     cvDestroyWindow("dst");27     return 0;28 }29 30 int Otsu(IplImage* src)31 {32     int height = src->height;33     int width = src->width;34     long size = height * width;35 36     //histogram    37     float histogram[256] = { 0 };38     for (int m = 0; m < height; m++)39     {40         unsigned char* p = (unsigned char*)src->imageData + src->widthStep * m;41         for (int n = 0; n < width; n++)42         {43             histogram[int(*p++)]++;44         }45     }46 47     int threshold;48     long sum0 = 0, sum1 = 0; //存储前景的灰度总和和背景灰度总和  49     long cnt0 = 0, cnt1 = 0; //前景的总个数和背景的总个数  50     double w0 = 0, w1 = 0; //前景和背景所占整幅图像的比例  51     double u0 = 0, u1 = 0;  //前景和背景的平均灰度  52     double variance = 0; //最大类间方差  53     int i, j;54     double u = 0;55     double maxVariance = 0;56     for (i = 1; i < 256; i++) //一次遍历每个像素  57     {58         sum0 = 0;59         sum1 = 0;60         cnt0 = 0;61         cnt1 = 0;62         w0 = 0;63         w1 = 0;64         for (j = 0; j < i; j++)65         {66             cnt0 += histogram[j];67             sum0 += j * histogram[j];68         }69 70         u0 = (double)sum0 / cnt0;71         w0 = (double)cnt0 / size;72 73         for (j = i; j <= 255; j++)74         {75             cnt1 += histogram[j];76             sum1 += j * histogram[j];77         }78 79         u1 = (double)sum1 / cnt1;80         w1 = 1 - w0; // (double)cnt1 / size;  81 82         u = u0 * w0 + u1 * w1; //图像的平均灰度  83         printf("u = %f\n", u);84         //variance =  w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);  85         variance = w0 * w1 *  (u0 - u1) * (u0 - u1);86         if (variance > maxVariance)87         {88             maxVariance = variance;89             threshold = i;90         }91     }92 93     printf("threshold = %d\n", threshold);94     return threshold;95 }
        
       
      
     

 

第二份代码：

1 #include 
     
        2 #include 
      
         3   4 int otsu(IplImage *image)  5 {  6     assert(NULL != image);  7   8     int width = image->width;  9     int height = image->height; 10     int x = 0, y = 0; 11     int pixelCount[256]; 12     float pixelPro[256]; 13     int i, j, pixelSum = width * height, threshold = 0; 14  15     uchar* data = (uchar*)image->imageData; 16  17     //初始化 18     for (i = 0; i < 256; i++) 19     { 20         pixelCount[i] = 0; 21         pixelPro[i] = 0; 22     } 23  24     //统计灰度级中每个像素在整幅图像中的个数 25     for (i = y; i < height; i++) 26     { 27         for (j = x; j 
       
        widthStep + j]]++; 30         } 31     } 32  33     //计算每个像素在整幅图像中的比例 34     for (i = 0; i < 256; i++) 35     { 36         pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum); 37     } 38  39     //经典ostu算法,得到前景和背景的分割 40     //遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值 41     float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0; 42     for (i = 0; i < 256; i++) 43     { 44         w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0; 45  46         for (j = 0; j < 256; j++) 47         { 48             if (j <= i) //背景部分 49             { 50                 //以i为阈值分类，第一类总的概率 51                 w0 += pixelPro[j]; 52                 u0tmp += j * pixelPro[j]; 53             } 54             else       //前景部分 55             { 56                 //以i为阈值分类，第二类总的概率 57                 w1 += pixelPro[j]; 58                 u1tmp += j * pixelPro[j]; 59             } 60         } 61  62         u0 = u0tmp / w0;        //第一类的平均灰度 63         u1 = u1tmp / w1;        //第二类的平均灰度 64         u = u0tmp + u1tmp;        //整幅图像的平均灰度 65         //计算类间方差 66         deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u); 67         //找出最大类间方差以及对应的阈值 68         if (deltaTmp > deltaMax) 69         { 70             deltaMax = deltaTmp; 71             threshold = i; 72         } 73     } 74     //返回最佳阈值; 75     return threshold; 76 } 77  78 int main(int argc, char* argv[]) 79 { 80     IplImage* srcImage = cvLoadImage("house.jpg", 0); 81     assert(NULL != srcImage); 82  83     cvNamedWindow("src"); 84     cvShowImage("src", srcImage); 85  86     IplImage* biImage = cvCreateImage(cvGetSize(srcImage), 8, 1); 87  88     //计算最佳阈值 89     int threshold = otsu(srcImage); 90     printf("threshold = %d\n", threshold); 91     //对图像二值化 92     cvThreshold(srcImage, biImage, threshold, 255, CV_THRESH_BINARY); 93  94     cvNamedWindow("binary"); 95     cvShowImage("binary", biImage); 96  97     cvWaitKey(0); 98  99     cvReleaseImage(&srcImage);100     cvReleaseImage(&biImage);101     cvDestroyWindow("src");102     cvDestroyWindow("binary");103 104     return 0;105 }
       
      
     

 

三、代码运行结果

代码1运行结果：

代码2运行结果：