OTSU算法
OTSU算法以最佳门限将图像灰度直方图分割成两部分,使两部分类间方差取最大值,即分离性最大。
设图像灰度级 ,第 级象素 个,总象素,则第级灰度出现的概率为 。
设灰度门限值为 ,则图像像素按灰度级被分为两类:
,
图像总平均灰度级:
类的平均灰度级为: ,像素数为:
类的平均灰度级为: , 像素数为:
两部分图像所占比例分别为:
对,均值作处理:
图像总均值可化为:
类间方差:
化为:
从变化,使 最大的 即为所求之最佳门限。
称为目标选择函数。
My code:
int i; //图像灰度级 from 0--255
double n[256]; //第i级像素n[i]个
double Pro[256]; //第i级灰度出现的概率
double Tavegray = 0; //图像总平均灰度级
double avegray0[256]; //C0类的平均灰度级
double avegray1[256]; //C1类的平均灰度级
int N = 0; //总像素
int N0 = 0; //C0类的像素数
int N1 = 0; //C1类的像素数
double w0 = 0, w1 =0; //分别为两部分图像所占比例
double mean0 = 0, mean1 =0, mean = 0; //对C0,C1 均值作处理,及图像总均值
double var[256]; //类间方差
int k = 0; //灰度门限值
double w[256]; //w[k]当门限值为k时,C0类图像所占的比值
int Threshold; //门限值
//数组初始化
for( i = 0; i <256; i++ )
{
n[i] = 0;
Pro[i] = 0;
avegray0[i] = 0;
avegray1[i] = 0;
var[i] = 0;
w[i] = 0;
}
int I;
LPSTRlpDibHdr = (LPSTR)::GlobalLock(hDib);
LPBITMAPINFOHEADERlpbmi = (LPBITMAPINFOHEADER)lpDibHdr; //得到指向位图的信息头结构的指针
LPSTRlpDibBits = lpDibHdr - lpbmi->biSize; //得到指向位图数据的指针
for( int y = 0; y <lpbmi->biHeight; y++ )
{
for( int x = 0; x < lpbmi->biWidth; x++ )
{
COLORREFc;
BYTE*pDate = (BYTE*)lpDibBits + (lpbmi->biHeight-y) * getpiont( hDib )+ ( x ) *3; //((lpbmi->biWidth*24+31)/32*4)
c= RGB( pDate[2], pDate[1], pDate[0]);
// COLORREF cNew;
// BYTE* pDateNew = (BYTE*)lpNewDibBits +(lpNewbmi->biHeight -y) * getpiont(hDib) + ( x ) * 3; //((lpbmi->biWidth*24+31)/32*4)
// cNew = RGB( pDateNew[2], pDateNew[1], pDateNew[0]);
I= ( pDate[2] + pDate[1] + pDate[0] ) / 3;
n[I]++;
}
}
//总象素
for( i = 0; i <256; i++ )
{
N += n[i];
}
//第i级灰度出现的概率
for( i = 0; i <256; i++ )
{
Pro[i] = ( n[i] / N );
}
//图像总平均灰度级
for( i = 1; i <=256; i++ )
{
Tavegray += ( i * Pro[i - 1] );
}
for ( k = 0; k< 256; k++ )
{
for( i = 1; i < k+1; i++ )
{
avegray0[k] += i*Pro[i];//C0类的平均灰度级
// N0 +=n[i]; //C0类像素数
w0 += Pro[i];
}
w[k]= w0;//C0部分图像所占比例
// avegray1[k]= Tavegray - avegray0[k];//C1类的平均灰度级和像素数
// N1= N - N0; //C1类像素数
w1= 1 - w[k]; //C0部分图像所占比例
if ( w[k] != 0 ) //对C0均值作处理
{
mean0 = avegray0[k]/ w[k];
}
else mean0 = 0;
mean1= ( Tavegray - avegray0[k] ) / ( 1 - w[k] );//对C1均值作处理
mean= ( w0 * mean0 ) + ( w1 * mean1 );
var[k]= (( ( mean * w[k] ) - avegray0[k]) * ( ( mean * w[k] ) - avegray0[k])) / (w[k] * ( 1 - w[k]));
}
int temp = 0;
for( k = 0; k< 256; k++ )
{
if ( var[k] >temp)
{
temp = var[k];
Threshold = k;
}
}
OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。下面的代码最早由Ryan Dibble提供,此后经过多人Joerg.Schulenburg,R.Z.Liu 等修改,补正。
算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。划分点就是求得的阈值。
parameter: *image --- buffer for image
rows, cols --- size of image
x0, y0, dx, dy --- region of vector used for computing threshold
vvv --- debug option, is 0, no debug information outputed
*/
/*======================================================================*/
/* OTSU global thresholdingroutine */
/* takes a 2D unsigned char array pointer, number of rows,and */
/* number of cols in the array. returns the value of thethreshold */
/*======================================================================*/
int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, intdy, int vvv)
{
unsigned char*np; // 图像指针
int thresholdValue=1; // 阈值
intihist[256]; // 图像直方图,256个点
int i, j,k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;
// 对直方图置零...
memset(ihist, 0, sizeof(ihist));
gmin=255; gmax=0;
// 生成直方图
for (i = y0 + 1; i < y0 + dy - 1; i++) {
np = &image[i*cols+x0+1];
for (j = x0 + 1; j < x0 + dx - 1; j++) {
ihist[*np]++;
if(*np > gmax) gmax=*np;
if(*np < gmin) gmin=*np;
np++; /* next pixel */
}
}
// set up everything
sum = csum = 0.0;
n = 0;
for (k = 0; k <= 255; k++) {
sum += (double) k * (double) ihist[k]; /*x*f(x) 质量矩*/
n +=ihist[k]; /* f(x) 质量 */
}
if (!n) {
// if n has no value, there is problems...
fprintf (stderr, "NOT NORMAL thresholdValue =160\n");
return (160);
}
// do the otsu global thresholdingmethod
fmax = -1.0;
n1 = 0;
for (k = 0; k < 255; k++) {
n1 += ihist[k];
if (!n1) { continue; }
n2 = n - n1;
if (n2 == 0) { break; }
csum += (double) k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax) {
fmax = sb;
thresholdValue = k;
}
}
// at this point we have ourthresholding value
// debug code to displaythresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%dgmax=%d\n",
thresholdValue, gmin, gmax);
return(thresholdValue);
}
matlab的OTSU算法;
function threshold=ostu(filename);
x=imread('filename');
%figure;
%imshow(x);
[m,n]=size(x);
N=m*n;
num=zeros(1,256);
p=zeros(1,256);
for i=1:m
for j=1:n
num(x(i,j)+1)=num(x(i,j)+1)+1;
end
end
for i=0:255;
p(i+1)=num(i+1)/N;
end
totalmean=0;
for i=0:255;
totalmean=totalmean+i*p(i+1);
end
maxvar=0;
for k=0:255
kk=k+1;
zerosth=sum(p(1:kk));
firsth=0;
for h=0:k
firsth=firsth+h*p(h+1);
end
var=totalmean*zerosth-firsth;
var=var*var;
var=var/(zerosth*(1-zerosth)+0.01);
var=sqrt(var);
if(var>maxvar)
maxvar=var;
point=k;
end
end
threshold=point;
网址:http://wenku.baidu.com/link?url=tOgi4aWIUlfxOSYWzXFRQnicPsTGDa048Yj9g_jySVmInlwnyhdCJp6qFVCQtYAqscDnnu7G4PCu-NDrw8uJ08Kt34gEptNX3GmcCWE-gE7