大津法由大津于1979年提出，对图像Image，记t为前景与背景的分割阈值，前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1。图像的总平均灰度为：u=w0*u0+w1*u1。从最小灰度值到最大灰度值遍历t，当t使得值g=w0*(u0-u)2+w1*(u1-u)2 最大时t即为分割的最佳阈值。对大津法可作如下理解：该式实际上就是类间方差值，阈值t分割出的前景和背景两部分构成了整幅图像，而前景取值u0，概率为 w0，背景取值u1，概率为w1，总均值为u，根据方差的定义即得该式。因方差是灰度分布均匀性的一种度量,方差值越大,说明构成图像的两部分差别越大, 当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小，因此使类间方差最大的分割意味着错分概率最小。
直接应用大津法计算量较大，因此我们在实现时采用了等价的公式g=w0*w1*(u0-u1)2。部分计算过程如下：
//遍历所有灰度值求Max g。
for intCurrentLevel:=0 to intArrLen do
begin
if intSclGrayLevel[intCurrentLevel]=0 then
continue
else
begin
//计算当阈值为intCurrentLevel时的g
intCount:=0;
intSumPels:=0;
for intLoop:=0 to intCurrentLevel do
begin
intCount:=intCount+intSclGrayLevel[intLoop];
intSumPels:=intSumPels+intSumPelsArr[intLoop];
end;
w0:=intCount/intSize;
u0:=intSumPels/intCount;
w1:=1-w0;
if intSize-intCount<>0 then
u1:=(intTotalPels-intSumPels)/(intSize-intCount)
else
u1:=0;
RlTempO:=w0*w1*(u0-u1)*(u0-u1);
if RlTempO>RlMaxO then
begin
RlMaxO:=RlTempO;
Result:=intCurrentLevel;
end;
end;
我们在测试中发现：大津法选取出来的阈值非常理想，对各种情况的表现都较为良好。虽然它在很多情况下都不是最佳的分割，但分割质量通常都有一定的保障，可以说是最稳定的分割。由上可知，大津算法是一种较为通用的分割算法。在它的思想的启迪下，人们进一步提出了多种类似的评估阈值的算法，具体可参加【5】、【6】等。
OTSU的算法，很好用，好不容易才找到的。
/*
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 thresholding routine */
/* takes a 2D unsigned char array pointer, number of rows, and */
/* number of cols in the array. returns the value of the threshold */
/*======================================================================*/
int otsu (unsigned char *image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
{
unsigned char *np; // 图像指针
int thresholdValue=1; // 阈值
int ihist[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 = ℑ[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 thresholding method
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 our thresholding value
// debug code to display thresholding values
if ( vvv & 1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d/n",
thresholdValue, gmin, gmax);
return(thresholdValue);
}
///
Otsu算法(大律法或最大类间方差法)
一、Otsu最大类间方差法原理
利用阈值将原图像分成前景，背景两个图象。
前景：用n1,csum,m1来表示在当前阈值下的前景的点数，质量矩，平均灰度
后景：用n2,sum-csum,m2来表示在当前阈值下的背景的点数，质量矩，平均灰度
当取最佳阈值时，背景应该与前景差别最大，关键在于如何选择衡量差别的标准，而在otsu算法中这个衡量差别的标准就是最大类间
方差（英文简称otsu，这也就是这个算法名字的来源），在本程序中类间方差用sb表示，最大类间方差用fmax
关于最大类间方差法（otsu）的性能:
类间方差法对噪音和目标大小十分敏感，它仅对类间方差为单峰的图像产生较好的分割效果。
当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好，但是类间方差法是用时最少的。
最大类间方差法（otsu）的公式推导：
记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的公式也就是下面程序中的sb的表达式。当方差g最大时，可以认为此时前景和背景差异最大，此时的灰
度t是最佳阈值sb= w1*w2*(u1-u0)*(u0-u1)
算法实现1：
unsafepublicintGetThreshValue(Bitmapimage)
{
BitmapDatabd=image.LockBits(new Rectangle(0,0,image.Width, image.Height), ImageLockMode.WriteOnly, image.PixelFormat);
byte*pt=(byte*)bd.Scan0;
int[] pixelNum=newint[256]; //图象直方图，共256个点
bytecolor;
byte*pline;
int n,n1,n2;
int total; //total 为总和，累计值
doublem1,m2,sum,csum,fmax, sb; //sb 为类间方差，fmax存储最大方差值
int k,t, q;
int threshValue =1; // 阈值
int step=1;
switch(image.PixelFormat)
{
casePixelFormat.Format24bppRgb:
step=3;
break;
casePixelFormat.Format32bppArgb:
step=4;
break;
casePixelFormat.Format8bppIndexed:
step=1;
break;
}
//生成直方图
for(inti =0;i
{
pline=pt+i *bd.Stride;
for(intj =0;j
{
color=*(pline+j *step); //返回各个点的颜色，以RGB表示
pixelNum[color]++; //相应的直方图加1
}
}
//直方图平滑化
for(k =0;k<=255;k++)
{
total =0;
for(t =-2;t <=2;t++) //与附近2个灰度做平滑化，t值应取较小的值
{
q=k+t;
if (q<0) //越界处理
q=0;
if (q>255)
q=255;
total =total +pixelNum[q]; //total 为总和，累计值
}
//平滑化，左边2个+中间1个+右边2个灰度，共5个，所以总和除以5，后面加0.5是用修正值
pixelNum[k]=(int)((float)total/ 5.0+0.5);
}
//求阈值
sum=csum=0.0;
n=0;
//计算总的图象的点数和质量矩，为后面的计算做准备
for(k =0;k<=255;k++)
{
//x*f(x)质量矩，也就是每个灰度的值乘以其点数（归一化后为概率），sum为其总和
sum+=(double)k *(double)pixelNum[k];
n+=pixelNum[k]; //n 为图象总的点数，归一化后就是累积概率
}
fmax=-1.0; //类间方差sb不可能为负，所以fmax初始值为-1不影响计算的进行
n1=0;
for(k =0;k<255;k++) //对每个灰度（从0到255）计算一次分割后的类间方差sb
{
n1+=pixelNum[k]; //n1 为在当前阈值遍前景图象的点数
if (n1==0){continue;} //没有分出前景后景
n2=n-n1; //n2 为背景图象的点数
//n2 为0表示全部都是后景图象，与n1=0情况类似，之后的遍历不可能使前景点数增加，所以此时可以退出循环
if (n2==0){break;}
csum+=(double)k*pixelNum[k]; //前景的“灰度的值*其点数”的总和
m1=csum/ n1; //m1 为前景的平均灰度
m2=(sum -csum)/ n2; //m2 为背景的平均灰度
sb=(double)n1*(double)n2*(m1-m2)*(m1- m2); //sb为类间方差
if (sb>fmax) //如果算出的类间方差大于前一次算出的类间方差
{
fmax=sb; //fmax 始终为最大类间方差（otsu）
threshValue =k; //取最大类间方差时对应的灰度的k就是最佳阈值
}
}
image.UnlockBits(bd);
image.Dispose();
returnthreshValue;
}
算法实现2：
Otsu算法步骤如下：
设图象包含L个灰度级(0,1…,L-1)，灰度值为i的的象素点数为Ni，图象总的象素点数为
N=N0+N1+...+N(L-1)。灰度值为i的点的概率为：P(i)= N(i)/N.
门限t将整幅图象分为暗区c1和亮区c2两类，则类间方差σ是t的函数：σ=a1*a2(u1-u2)^2(2)
式中，aj为类cj的面积与图象总面积之比，a1=sum(P(i))i->t,a2= 1-a1;
uj为类cj的均值，u1= sum(i*P(i))/a10->t,u2= sum(i*P(i))/a2,t+1->L-1,该法选择最佳门限t^
使类间方差最大，即：令Δu=u1-u2，σb = max{a1(t)*a2(t)Δu^2}
代码实现：
int otsu(IplImage*image,int rows,int cols,int x0,int y0,int dx,int dy,int vvv)
{
unsignedchar*np;// 图像指针
int thresholdValue=1;// 阈值
int ihist[256]; // 图像直方图，256个点
int i, j,k;// variouscounters
int n,n1,n2,gmin,gmax;
doublem1,m2,sum,csum,fmax,sb;
// 对直方图置零
memset(ihist,0,sizeof(ihist));
gmin=255;gmax=0;
// 生成直方图
for(j=y0;j
{
for(i=0;i
{
unsignedchartemp=CV_IMAGE_ELEM(image,uchar,j,i);
ihist[temp]++;
}
}
//setupeverything
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 nhasnovalue,thereis problems
fprintf(stderr,"NOTNORMALthresholdValue=160\n");
return(160);
}
//dotheotsuglobalthresholdingmethod
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);
if (sb>fmax){
fmax=sb;
thresholdValue=k;
}
}
//atthis pointwehaveourthresholdingvalue，debugcodetodisplaythresholdingvalues
if (vvv&1)
fprintf(stderr,"#OTSU:thresholdValue=%dgmin=%dgmax=%d\n",thresholdValue,gmin,gmax);
return(thresholdValue);
}
简称OTSU。它是按图像的灰度特性,将图像分成背景和目标2部分。背景和目标之间的类间方差
越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部
分差别变小。因此,使类间方差最大的分割意味着错分概率最小。
对于图像I(x,y),前景(即目标)和背景的分割阈值记作T,属于前景的像素点数占整幅图像的比
例记为ω0,其平均灰度μ0;背景像素点数占整幅图像的比例为ω1,其平均灰度为μ1。图像的总平均
灰度记为μ,类间方差记为g。
假设图像的背景较暗,并且图像的大小为M×N,
图像中像素的灰度值小于阈值T的像素个数记作N0,像素灰度大于阈值T的像素个数记作N1,则有:
ω0=N0/M×N
2100
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
