记t为前景与背景的分割阈值,前景点数占图像比例为w0,平均灰度为u0;背景点数占图像比例为w1,平均灰度为u1。
则图像的总平均灰度为:u=w0*u0+w1*u1。
前景和背景图象的方差:g=w0*(u0-u)(u0-u)+w1(u1-u)(u1-u)=w0*w1(u0-u1)*(u0-u1),此公式为方差公式。
当方差g最大时,可以认为此时前景和背景差异最大,此时的灰度t是最佳阈值sb = w0*w1*(u1-u0)*(u0-u1)
//===============================================================================
//
// 函数名称:OTSU_threshold
// 功能说明:经典大津算法 动态阈值
// 修改时间:2016-8-31
// 备 注:
// Otsu实现思路
// 1.计算0~255各灰阶对应的像素个数,保存至一个数组中,该数组下标是灰度值,保存内容是当前灰度值对应像素数
// 2.计算背景图像的平均灰度、背景图像像素数所占比例
// 3.计算前景图像的平均灰度、前景图像像素数所占比例
// 4.遍历0~255各灰阶,计算并寻找类间方差极大值,方差是实际值与期望值之差平方的平均值
//***ps.一幅图像所谓的前景背景指的是在遍历所以灰度的过程中,以当前的灰度为界,分前景,背景***
//===============================================================================
#define Gourd 256
uint8 OTSU_threshold(uint8 *pic,uint16 num)
{
uint16 i=0;
uint16 Histogram[Gourd];//直方图histogram
for (i=0;i<Gourd;i++)
Histogram[i]=0;//数组清零
for (i=0;i<num;i++)
{
Histogram[(int)pic[i]*Gourd/256]++;//遍历每个像素,计算每个灰度级的像素个数和。
}
float pt[Gourd],w[Gourd],u[Gourd],o[Gourd],Ut;
pt[0]=(float)Histogram[0]/num;
w[0]=pt[0];
u[0]=w[0];
for(i=1;i<Gourd;i++)
{
pt[i]=(float)Histogram[i]/num; //灰度级为i的像素个数占总像素的比例
w[i]=w[i-1]+pt[i];//进行pt[i]的累加。对图像像素数所占比例进行累加
u[i]=u[i-1]+i*pt[i];//进行i*pt[i]的累加,从而计算出平均灰度(每个灰度级*每个灰度级所占的比例)
};
Ut=u[Gourd-1];//整幅图像平均灰度
for(i=0;i<Gourd;i++)
{
o[i]=(1-pt[i])*(u[i]*u[i]/w[i]+(u[i]-Ut)*(u[i]-Ut)/(1-w[i]));//方差
};
int maxi=0;
float maxo=0;
for(i=0;i<Gourd;i++)
{
if(o[i]!=0x7FC0000)
if(o[i]>maxo){maxo=o[i];maxi=i;}//遍历0~255各灰阶,计算并寻找类间方差极大值,当找到时,i为对应类间方差的灰度值
}
return maxi*256/Gourd;
}
//Otsu的C++实现
//转自
> http://www.cnblogs.com/gzy-zju-edu/articles/4202252.html
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;
}