正文
最近在学习图像分割反面的知识,在冈萨雷斯的那本书上看到Otsu算法,身边的同学都是用Matlab来实现这个算法。我觉得Matlab写得话,但是代码的效率应该不会高。于是又恶补了一些OpenCV的一些基本知识,然后看了Augusdi的博客,分析了一下他的代码并附上,第二份代码是来自某一位大牛的代码,写得更加易懂。
算法的介绍
otsu法(最大类间方差法,有时也称之为大津算法)使用的是聚类的思想,把图像的灰度数按灰度级分成2个部分,使得两个部分之间的灰度值差异最大,每个部分之间的灰度差异最小,通过方差的计算来寻找一个合适的灰度级别 来划分。 所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法,计算简单,不受图像亮度和对比度的影响。因此,使类间方差最大的分割意味着错分概率最小。
设t为设定的阈值。
w0
: 分开后 前景像素点数占图像的比例
u0
: 分开后 前景像素点的平均灰度
w1
:分开后 被景像素点数占图像的比例
u1
: 分开后 被景像素点的平均灰度
u=w0∗u0+w1∗u1
:图像总平均灰度
从L个灰度级遍历t,使得t为某个值的时候,前景和背景的方差最大, 则 这个 t 值便是我们要求得的阈值。其中,方差的计算公式如下:
此公式计算量较大,可以采用:
由于Otsu算法是对图像的灰度级进行聚类,so 在执行otsu算法之前,需要计算该图像的灰度直方图。
代码实现算法
第一份代码:来自Augusdi
#include<iostream>
#include<opencv\cv.h>
#include<opencv2\highgui.hpp>
int Otsu(IplImage* src);
int main()
{
IplImage* img = cvLoadImage("E:\\pp.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;
}
第二份代码:
#include <opencv\cv.h>
#include <opencv2\highgui\highgui.hpp>
int otsu(IplImage *image)
{
assert(NULL != image);
int width = image->width;
int height = image->height;
int x = 0, y = 0;
int pixelCount[256];
float pixelPro[256];
int i, j, pixelSum = width * height, threshold = 0;
uchar* data = (uchar*)image->imageData;
//初始化
for (i = 0; i < 256; i++)
{
pixelCount[i] = 0;
pixelPro[i] = 0;
}
//统计灰度级中每个像素在整幅图像中的个数
for (i = y; i < height; i++)
{
for (j = x; j <width; j++)
{
pixelCount[data[i * image->widthStep + j]]++;
}
}
//计算每个像素在整幅图像中的比例
for (i = 0; i < 256; i++)
{
pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum);
}
//经典ostu算法,得到前景和背景的分割
//遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值
float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;
for (i = 0; i < 256; i++)
{
w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;
for (j = 0; j < 256; j++)
{
if (j <= i) //背景部分
{
//以i为阈值分类,第一类总的概率
w0 += pixelPro[j];
u0tmp += j * pixelPro[j];
}
else //前景部分
{
//以i为阈值分类,第二类总的概率
w1 += pixelPro[j];
u1tmp += j * pixelPro[j];
}
}
u0 = u0tmp / w0; //第一类的平均灰度
u1 = u1tmp / w1; //第二类的平均灰度
u = u0tmp + u1tmp; //整幅图像的平均灰度
//计算类间方差
deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);
//找出最大类间方差以及对应的阈值
if (deltaTmp > deltaMax)
{
deltaMax = deltaTmp;
threshold = i;
}
}
//返回最佳阈值;
return threshold;
}
int main(int argc, char* argv[])
{
IplImage* srcImage = cvLoadImage("E:\\pp.jpg", 0);
assert(NULL != srcImage);
cvNamedWindow("src");
cvShowImage("src", srcImage);
IplImage* biImage = cvCreateImage(cvGetSize(srcImage), 8, 1);
//计算最佳阈值
int threshold = otsu(srcImage);
//对图像二值化
cvThreshold(srcImage, biImage, threshold, 255, CV_THRESH_BINARY);
cvNamedWindow("binary");
cvShowImage("binary", biImage);
cvWaitKey(0);
cvReleaseImage(&srcImage);
cvReleaseImage(&biImage);
cvDestroyWindow("src");
cvDestroyWindow("binary");
return 0;
}
代码运行结果