图像处理之基于OTSU阈值二值化
这篇文章主要参考了:http://blog.****.net/liyuanbhu/article/details/49387483 还有http://blog.****.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 给出的类间方差定义为:
那么这个最佳的阈值t 就是使得 ICV 最大的那个值。
对于上面的测试图像,我们可以遍历 t 的各种取值,计算 ICV。之后可以画出这样的ICV 曲线(绿色线条):
在大津法中,我们定义组内方差为
通过选择使得上述组内方差最小化时的阈值 t,就可以使得图像中的前景和背景尽可能的被区别开(假设我们将最终图像里被分开的两部分称为前景和背景)。w0和w1分别是一个像素可能属于前景或背景的概率,而 σ 表示两个类别的方差。如果一个图像的直方图有L个等级(一般L=256),那么在给定阈值 t的情况下,w0和w1分别定义为
大津展之证明最小化组内方差(intra-class variance)与最大化组间方差(inter-class variance)是等价的,于是有
又因为(其中 μ 表示均值或期望)
可以推出
这个证明仅仅涉及一些算术上的推导,我简单演示如下
- #include<opencv2\opencv.hpp>
- #include<cmath>
- #include<iostream>
- using namespace cv;
- using namespace std;
- //OTSU算法函数实现
- int OTSU(Mat &srcImage)
- {
- int nRows = srcImage.rows;
- int nCols = srcImage.cols;
- int threshold = 0;
- double max = 0.0;
- double AvePix[256];
- int nSumPix[256];
- double nProDis[256];
- double nSumProDis[256];
- for (int i = 0; i < 256; i++)
- {
- AvePix[i] = 0.0;
- nSumPix[i] = 0;
- nProDis[i] = 0.0;
- nSumProDis[i] = 0.0;
- }
- for (int i = 0; i < nRows; i++)
- {
- for (int j = 0; j < nCols; j++)
- {
- nSumPix[(int)srcImage.at<uchar>(i, j)]++;
- }
- }
- for (int i = 0; i < 256; i++)
- {
- nProDis[i] = (double)nSumPix[i] / (nRows*nCols);
- }
- AvePix[0] = 0;
- nSumProDis[0] = nProDis[0];
- for (int i = 1; i < 256; i++)
- {
- nSumProDis[i] = nSumProDis[i - 1] + nProDis[i];
- AvePix[i] = AvePix[i - 1] + i*nProDis[i];
- }
- double mean = AvePix[255];
- for (int k = 1; k < 256; k++)
- {
- double PA = nSumProDis[k];
- double PB = 1 - nSumProDis[k];
- double value = 0.0;
- if (fabs(PA) > 0.001 && fabs(PB) > 0.001)
- {
- double MA = AvePix[k];
- double MB = (mean - PA*MA)/PB;
- value = (double)(PA * PB * pow((MA - MB), 2));
- //或者这样value = (double)(PA * PB * pow((MA-MB),2));//类间方差
- //pow(PA,1)* pow((MA - mean),2) + pow(PB,1)* pow((MB - mean),2)
- if (value > max)
- {
- max = value;
- threshold = k;
- }
- }
- }
- return threshold;
- }
- int main()
- {
- Mat srcImage = imread("img.jpg");
- if (!srcImage.data)
- {
- printf("could not load image...\n");
- return -1;
- }
- Mat srcGray;
- cvtColor(srcImage, srcGray, CV_BGR2GRAY);
- imshow("srcGray", srcGray);
- //调用二值化函数得到最佳阈值
- int otsuThreshold = OTSU(srcGray);
- cout << otsuThreshold << endl;//输出最佳阈值
- Mat otsuResultImage = Mat::zeros(srcGray.rows, srcGray.cols, CV_8UC1);
- //利用得到的阈值进行二值操作
- for (int i = 0; i < srcGray.rows; i++)
- {
- for (int j = 0; j < srcGray.cols; j++)
- {
- if (srcGray.at<uchar>(i, j) > otsuThreshold)
- {
- otsuResultImage.at<uchar>(i, j) = 255;
- }
- else
- {
- otsuResultImage.at<uchar>(i, j) = 0;
- }
- }
- }
- imshow("otsuResultImage", otsuResultImage);
- waitKey(0);
- return 0;
- }
效果图: