基本思想
OTSU算法主要用于在灰度图像中查找一个阈值k,该阈值把图像区分为前景图像C0和背景图像C1。
关键在于阈值是怎样确定的。OTSU算法的思路类似于聚类,聚类就是给定一定的数据,让你根据数据的特性进行分类,使得这两类的数据能够分类合理,即常说的无监督学习。相对应的监督学习则是给定数据及分类对训练器进行训练然后对新数据进行估计。
对于这个阈值k,我们可以得到两个分类对应的概率 w0, w1 和 均值 u0, u1。
以及对应的类方差。
而在分类的过程中,通常我们需要采用一个标准来评判阈值是否合适。经过OTSU在论文中的介绍,采用的是最大类间方差。
具体可以查看论文《[1979 IEEE] OTSU A threshold selection method from gray-level histograms》
最大类间方差:
代码实现
以下为获取阈值的代码:
<span style="font-size:14px;font-weight: normal;">int getOtsuTh(unsigned char buff[],long size)
{
int i,j;
unsigned char indexMax=0;
double w0;
double u0;
double w1;
double u1;
double u;
unsigned int cnt0;
unsigned int cnt1;
unsigned int sum1;
unsigned int sum0;
double vaMax=-1;
double va;
jint* graylvl=(jint*)malloc(256*sizeof(jint));
memset(graylvl,0,256*sizeof(jint));
// 灰度值统计
for(i=0;i<size;i++)
{
graylvl[buff[i]]++;
}
for(i=1;i<256;i++)
{
cnt0=0;
cnt1=0;
sum0=0;
sum1=0;
for(j=0;j<i;j++)
{ // 对类C0进行灰度值统计
cnt0+=graylvl[j];
sum0+=graylvl[j]*j;
}
w0= (double)cnt0/size;
u0= (double)sum0/cnt0;
for(j=i;j<256;j++)
{ // 对类C1进行灰度值统计
cnt1+=graylvl[j];
sum1+=graylvl[j]*j;
}
w1= 1-w0;
u1= (double)sum1/cnt1;
u= w0*u0+w1*u1;
va = w0*w1*(u0-u1)*(u0-u1);
if(va>vaMax)
{
vaMax=va;
indexMax = i;
}
}
free(graylvl);
return indexMax;
}</span>
该代码主要参考:http://blog.csdn.net/timidsmile/article/details/8493468
整个工程代码:下载地址
代码中使用的是图片的bin格式,即文件保存的是图像的灰度值,可参考《【图像处理】matlab操作图像bin文件》
<span style="font-size:14px;font-weight: normal;">int getOtsuTh(unsigned char buff[],long size)
{
int i,j;
unsigned char indexMax=0;
double w0;
double u0;
double w1;
double u1;
double u;
unsigned int cnt0;
unsigned int cnt1;
unsigned int sum1;
unsigned int sum0;
double vaMax=-1;
double va;
jint* graylvl=(jint*)malloc(256*sizeof(jint));
memset(graylvl,0,256*sizeof(jint));
// 灰度值统计
for(i=0;i<size;i++)
{
graylvl[buff[i]]++;
}
for(i=1;i<256;i++)
{
cnt0=0;
cnt1=0;
sum0=0;
sum1=0;
for(j=0;j<i;j++)
{ // 对类C0进行灰度值统计
cnt0+=graylvl[j];
sum0+=graylvl[j]*j;
}
w0= (double)cnt0/size;
u0= (double)sum0/cnt0;
for(j=i;j<256;j++)
{ // 对类C1进行灰度值统计
cnt1+=graylvl[j];
sum1+=graylvl[j]*j;
}
w1= 1-w0;
u1= (double)sum1/cnt1;
u= w0*u0+w1*u1;
va = w0*w1*(u0-u1)*(u0-u1);
if(va>vaMax)
{
vaMax=va;
indexMax = i;
}
}
free(graylvl);
return indexMax;
}</span>
后续将再补充对该代码的优化。