最大类间方差法

最大类间方差法是由日本学者大津于1979年提出的,是一种自适应的阈值确定的方法,又叫大津
法,简称OTSU。

它是按图像的灰度特性,将图像分成背景和目标2部分。

背景和目标之间的类间方差
越大,说明构成图像的2部分的差别越大,当部分目标错分为背景或部分背景错分为目标都会导致2部
分差别变小。

因此,使类间方差最大的分割意味着错分概率最小。

对于图像I(x,y),前景(即目标)和背景的分割阈值记作T,属于前景的像素点数占整幅图像的比
例记为ω0,其平均灰度μ0;背景像素点数占整幅图像的比例为ω1,其平均灰度为μ1。

图像的总平均
灰度记为μ,类间方差记为g。

假设图像的背景较暗,并且图像的大小为M×N,
图像中像素的灰度值小于阈值T的像素个数记作N0,像素灰度大于阈值T的像素个数记作N1,则有: ω0=N0/ M×N (1)
ω1=N1/ M×N (2)
N0+N1=M×N (3)
ω0+ω1=1 (4)
μ=ω0*μ0+ω1*μ1 (5)
g=ω0(μ0-μ)^2+ω1(μ1-μ)^2 (6)
将式(5)代入式(6),得到等价公式: g=ω0ω1(μ0-μ1)^2 (7)
采用遍历的方法得到使类间方差最大的阈值T,即为所求。

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))/a1 0->t,
u2 = sum(i*P(i))/a2, t+1->L-1
该法选择最佳门限t^ 使类间方差最大，即：令Δu=u1-u2，σb = max{a1(t)*a2(t)Δu^2}
OTSU方法计算图像二值化的自适应阈值
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 = &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++;
}
}
// set up everything
sum = csum = 0.0;
n = 0;
for (k = 0; k <= 255; k++) {
sum += (double) k * (double) ihist[k];
n += ihist[k];
}
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); 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);
}
///////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
另外一个
///////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
unsafe public int GetThreshValue(Bitmap image)
{
BitmapData bd = image.LockBits(new Rectangle(0, 0, image.Width, image.Height), ImageLockMode.WriteOnly, image.PixelFormat);
byte* pt = (byte*)bd.Scan0;
int[] pixelNum = new int[256]; //图象直方图，共256个点
byte color;
byte* pline;
int n, n1, n2;
int total; //total为总和，累计值
double m1, m2, sum, csum, fmax, sb; //sb为类间方差，fmax存储最大方差值
int k, t, q;
int threshValue = 1; // 阈值
int step = 1;
switch (image.PixelFormat)
{
case PixelFormat.Format24bppRgb:
step = 3;
break;
case PixelFormat.Format32bppArgb:
step = 4;
break;
case PixelFormat.Format8bppIndexed:
step = 1;
break;
}
//生成直方图
for (int i = 0; i < image.Height; i++)
{
pline = pt + i * bd.Stride;
for (int j = 0; j < image.Width; 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为总和，累计值
}
pixelNum[k] = (int)((float)total / 5.0 + 0.5); //平滑化，左边2个+中间1个+右边2个灰度，共5个，所以总和除以5，后面加0.5是用修正值
}
//求阈值
sum = csum = 0.0;
n = 0;
//计算总的图象的点数和质量矩，为后面的计算做准备
for (k = 0; k <= 255; k++)
{
sum += (double)k * (double)pixelNum[k]; //x*f(x)质量矩，也就是每个灰度的值乘以其点数（归一化后为概率），sum为其总和
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为背景图象的点数
if (n2 == 0) { break; } //n2为0表示全部都是后景图象，与n1=0情况类似，之后的遍历不可能使前景点数增加，所以此时可以退出循环
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();
return threshValue;
}。