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 算法是对图像的灰度级进行聚类, so 在执行 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 _tmain(int argc, _TCHAR* argv[])
- {
- IplImage* img = cvLoadImage("c:\\aSa.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 height=src->height;
- int width=src->width;
- long size = height * width;
- //histogram
- float histogram[256] = {0};
- for(int m=0; m <height; m++)
- {
- unsigned char* p=(unsigned char*)src->imageData + src->widthStep * m;
- for(int n = 0; n <width; n++)
- {
- histogram[int(*p++)]++;
- }
- }
- 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;
- 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;
- }
把 w1 写成 w0 .. 害我 debug 了好久~~ 总是不认真, 脑袋浑浑噩噩的... 这都看不出来....
==================
对了, 之前搜集的一个 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 _tmain(int argc, _TCHAR* argv[]) {
- IplImage* img = cvLoadImage("c:\\aSa.jpg",0); IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1); int threshold = Otsu(img); printf("threshold = %d\n", threshold); 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 height=src->height; int width=src->width; //histogram float histogram[256] = {
- 0
- }; for(int i=0; i <height; i++) {
- unsigned char* p=(unsigned char*)src->imageData + src->widthStep * i; for(int j = 0; j <width; j++) {
- histogram[*p++]++;
- }
- } //normalize histogram int size = height * width; for(int i = 0; i < 256; i++) {
- histogram[i] = histogram[i] / size;
- } //average pixel value float avgValue=0; for(int i=0; i < 256; i++) {
- avgValue += i * histogram[i]; // 整幅图像的平均灰度
- } int threshold; float maxVariance=0; float w = 0, u = 0; for(int i = 0; i < 256; i++) {
- w += histogram[i]; // 假设当前灰度 i 为阈值, 0~i 灰度的像素 (假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例 u += i * histogram[i]; // 灰度 i 之前的像素(0~i) 的平均灰度值: 前景像素的平均灰度值 float t = avgValue * w - u; float variance = t * t / (w * (1 - w) ); if(variance> maxVariance) {
- maxVariance = variance; threshold = i;
- }
- } return threshold;
- }
来源: http://www.bubuko.com/infodetail-2981046.html