图像处理之基于OTSU阈值二值化

这篇文章主要参考了：http://blog.csdn.net/liyuanbhu/article/details/49387483 还有http://blog.csdn.net/jia20003/article/details/43915983

机器视觉领域许多算法都要求先对图像进行二值化。这种二值化操作阈值的选取非常重要。阈值选取的不合适，可能得到的结果就毫无用处。今天就来讲讲一种自动计算阈值的方法。这种方法被称之为Otsu法。发明人是个日本人，叫做Nobuyuki Otsu (大津展之）。 
简单的说，这种算法假设一副图像由前景色和背景色组成，通过统计学的方法来选取一个阈值，使得这个阈值可以将前景色和背景色尽可能的分开。或者更准确的说是在某种判据下最优。与数理统计领域的 fisher 线性判别算法其实是等价的。

otsu算法中这个判据就是最大类间方差 (intra-class variance or the variance within the class)。下面就来详细说说什么是 intra-class variance。

我们知道一副灰度图像，可以计算它的颜色平均值，或者更进一步。可以计算出灰度直方图

比如下面的例子图片： 

这个图片拍摄的是一个条形码。在这个图中，前景色就是黑色的条形码，背景色是其余部分的灰色。那么我们可以计算出这个图像的灰度直方图。

 
图中那个大的峰是背景色的部分，小的峰是前景色。

灰度值的均值是 122. 我们称这个均值为 M。

现在任意选取一个灰度值 t，则可以将这个直方图分成前后两部分。我们称这两部分分别为 A 和 B。对应的就是前景色和背景色。这两部分各自的平均值成为 MA 和 MB。 
A 部分里的像素数占总像素数的比例记作 PA，B部分里的像素数占总像素数的比例记作 PB。 
Nobuyuki Otsu 给出的类间方差定义为： 

ICV=PA∗(MA−M)2+PB∗(MB−M)2

那么这个最佳的阈值 t  就是使得  ICV  最大的那个值。 
对于上面的测试图像，我们可以遍历 t 的各种取值，计算 ICV。之后可以画出这样的ICV 曲线（绿色线条）： 

在大津法中，我们定义组内方差为

通过选择使得上述组内方差最小化时的阈值 t，就可以使得图像中的前景和背景尽可能的被区别开（假设我们将最终图像里被分开的两部分称为前景和背景）。w0和w1分别是一个像素可能属于前景或背景的概率，而 σ 表示两个类别的方差。如果一个图像的直方图有L个等级（一般L=256），那么在给定阈值 t的情况下，w0和w1分别定义为

大津展之证明最小化组内方差（intra-class variance）与最大化组间方差（inter-class variance）是等价的，于是有

又因为（其中 μ 表示均值或期望）

可以推出

这个证明仅仅涉及一些算术上的推导，我简单演示如下

[cpp]  view plain  copy

1. #include<opencv2\opencv.hpp>  
2. #include<cmath>  
3. #include<iostream>  
4. using namespace cv;  
5. using namespace std;  
6. //OTSU算法函数实现  
7. int OTSU(Mat &srcImage)  
8. {  
9.     int nRows = srcImage.rows;
10.     int nCols = srcImage.cols;
11.   
12.     int threshold = 0;
13.     double max = 0.0;
14.     double AvePix[256];
15.     int nSumPix[256];
16.     double nProDis[256]; 
17.     double nSumProDis[256];
18.   
19.   
20.     
21.     for (int i = 0; i < 256; i++)  
22.     {  
23.         AvePix[i] = 0.0;  
24.         nSumPix[i] = 0;  
25.         nProDis[i] = 0.0;  
26.         nSumProDis[i] = 0.0;  
27.     }  
28.   
29.     for (int i = 0; i < nRows; i++)  
30.     {  
31.         for (int j = 0; j < nCols; j++)  
32.         {  
33.             nSumPix[(int)srcImage.at<uchar>(i, j)]++;  
34.         }  
35.     }  
36.   
37.     for (int i = 0; i < 256; i++)  
38.     {  
39.         nProDis[i] = (double)nSumPix[i] / (nRows*nCols);  
40.   
41.     }  
42.   
43.     
44.     AvePix[0] = 0;  
45.     nSumProDis[0] = nProDis[0];  
46.   
47.     
48.     for (int i = 1; i < 256; i++)  
49.     {  
50.         nSumProDis[i] = nSumProDis[i - 1] + nProDis[i];
51.         AvePix[i] = AvePix[i - 1] + i*nProDis[i];
52.     }  
53.   
54.     double mean = AvePix[255];
55.   
56.   
57.     for (int k = 1; k < 256; k++)  
58.     {  
59.         double PA = nSumProDis[k];
60.         double PB = 1 - nSumProDis[k];
61.         double value = 0.0;
62.         if (fabs(PA) > 0.001 && fabs(PB) > 0.001)  
63.         {  
64.             double MA = AvePix[k]; 
65.             double MB = (mean - PA*MA)/PB;
66.             value = (double)(PA * PB * pow((MA - MB), 2));
67.             //或者这样value = (double)(PA * PB * pow((MA-MB),2));//类间方差  
68.             //pow(PA,1)* pow((MA - mean),2) + pow(PB,1)* pow((MB - mean),2)  
69.             if (value > max)  
70.             {  
71.                 max = value;
72.                 threshold = k;
73.             }  
74.         }  
75.     }  
76.     return threshold;  
77. }  
78. int main()  
79. {  
80.     Mat srcImage = imread("img.jpg");  
81.     if (!srcImage.data)  
82.     {  
83.         printf("could not load image...\n");  
84.         return -1;  
85.     }  
86.     Mat srcGray;  
87.     cvtColor(srcImage, srcGray, CV_BGR2GRAY);
88.     imshow("srcGray", srcGray);  
89.     //调用二值化函数得到最佳阈值  
90.     int otsuThreshold = OTSU(srcGray);  
91.     cout << otsuThreshold << endl;//输出最佳阈值  
92.   
93.     Mat otsuResultImage = Mat::zeros(srcGray.rows, srcGray.cols, CV_8UC1);
94.   
95.     //利用得到的阈值进行二值操作  
96.     for (int i = 0; i < srcGray.rows; i++)  
97.     {  
98.         for (int j = 0; j < srcGray.cols; j++)  
99.         {  
100.             if (srcGray.at<uchar>(i, j) > otsuThreshold)  
101.             {  
102.                 otsuResultImage.at<uchar>(i, j) = 255;  
103.             }  
104.             else  
105.             {  
106.                 otsuResultImage.at<uchar>(i, j) = 0;  
107.             }  
108.         }  
109.     }  
110.     imshow("otsuResultImage", otsuResultImage);  
111.     waitKey(0);  
112.     return 0;  
113. }  

原图：

效果图：