一维最大类间方差法

最大类间方差法是由日本学者大津于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 countersint 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 everythingsum = csum = 0。

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：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为总和，累计值}//平滑化，左边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();return threshValue;}算法实现2：Otsu算法步骤如下：设图象包含L个灰度级(0,1…,L-1)，灰度值为i的的象素点数为Ni ，图象总的象素点数为N=N0+N1+...+N(L-1)。

由Otsu(大津展之)于1978年提出的最大类间方差法，是引起较多关注的一种阈值选取方法。

它是在判决分析或最小二乘原理的基础上推导出来的。

参考文献：[1] Otsu N. A threshold selection method from gray-level histogram. IEEE Trans,1979;SMC-9;62-66 下载地址算法思想：假设一幅图像有L个灰度级[1,2,…,L]。

灰度级为i的像素点的个数为n i，那么总的像素点个数就应该为N=n1+n2+…+n L。

为了讨论方便，我们使用归一化的灰度级直方图并且视为这幅图像的概率分布:现在假设我们通过一个灰度级为k的门限将这些像素点划分为两类:C0和C1(背景和目标,或者反之亦然)；C0表示灰度级为[1,…,k]的像素点，C1表示灰度级为[k+1,…,L]的像素点。

那么，每一类出现的概率以及各类的平均灰度级分别由下面的式子给出：以及其中，分别为灰度级从1到k的累积出现概率和平均灰度级（一阶累积矩），而是整幅图像的平均灰度级。

我们可以很容易验证，对于任意选定的k，都有：这两类的类内方差由下面的公式给出：这需要二阶累积矩(second-order cumulative moment，统计学概念)。

为了评价(灰度级k)这个门限“好”的程度，我们需要引入判别式分析中使用的判别式标准来测量（类的分离性测量）：其中：又根据式(9)，可以得出：这三个式子分别是类内方差、类间方差和灰度级的总方差。

然后，我们的问题就简化为一个优化问题，即寻找一个门限k使(12)式中给出的一个目标函数取最大值。

这个观点是出于这样一个猜想，一个好的阈值将会把灰度级分为两类，那么反过来说，就是如果一个门限能够在灰度级上将图像分割为最好的两类的话，那么这个门限就是最好的门限。

上面给出的判别式标准是分别求取λ、κ和η的最大值。

然而，对于κ而言，它又等于另外一个，比如κ=λ+1；而对于λ而言，又有η=λ/(λ+1)，因为始终存在下面的基本关系：效果展示：算法评价：就最大类间方差算法而言，优点是算法简单，当目标与背景的面积相差不大时，能够有效地对图像进行分割。

阈值确定⽅法⼀、问题重述图形（或图像）在计算机⾥主要有两种存储和表⽰⽅法。

⽮量图是使⽤点、直线或多边形等基于数学⽅程的⼏何对象来描述图形，位图则使⽤像素来描述图像。

⼀般来说，照⽚等相对杂乱的图像使⽤位图格式较为合适，⽮量图则多⽤于⼯程制图、标志、字体等场合。

⽮量图可以任意放缩，图形不会有任何改变。

⽽位图⼀旦放⼤后会产⽣较为明显的模糊，线条也会出现锯齿边缘等现象。

⽮量图从本质上只是使⽤曲线⽅程对图形进⾏的精确描述，在以像素为基本显⽰单元的显⽰器或打印机上是⽆法直接表现的。

将⽮量图转换成以像素点阵来表⽰的信息，再加以显⽰或打印，这个过程称之为栅格化（Rasterization），见图1。

栅格化的逆过程相对⽐较困难。

假设有⼀个形状较为简单的图标，保存成⼀定分辨率的位图⽂件。

我们希望将其⽮量化，请你建⽴合理的数学模型，尽量准确地提取出图案的边界线条，并将其⽤⽅程表⽰出来。

⼆、问题分析本题的要求是完成位图的⽮量化，通过建⽴合理的数学模型，将⼀个有⼀定分辨率的位图⽂件尽量准确地提取出图案的边界线条，最终将位图⽤⽅程的形式表⽰出来。

解决本问题的流程图见下图。

⾸先，通过MATLAB读取位图的各个像素的像素值（0-1），得到位图各个点的灰度值，通过最⼤类间⽅差法和最⼤熵法确定阈值，完成灰度的⼆值化，使各个像素点的灰度值全部由0或1表⽰。

其次，将位图的轮廓通过合适的算法提取出来，根据特征值对轮廓进⾏拟合。

最后，根据拟合的函数完成位图的⽮量图，完成其⽮量化过程，并通过对⽐⽮量图和原始位图对应的。

三、问题假设及符号说明3.1问题假设3.2符号说明四、模型建⽴4.1模型准备本题要求将⼀个形状较为简单的图标，保存成⼀定分辨率的位图⽂件，即将位图⽮量化。

阈值：指释放⼀个⾏为反应所需要的最⼩刺激强度，本⽂指像素点灰度值⼆值化的临界值。

4.2阈值的确定⽅法 4.2.1最⼤类间⽅差法最⼤类间⽅差法的基本思想是将待分割图像看作是由两类组成的整体，⼀类是背景，⼀类是⽬标[6]。

Otsu阈值最大类间方差法是一种确定阈值的算法，通过将图像分成前景和背景两部分，以最大化这两部分之间的方差。

这种方法是在判决分析最小二乘法原理的基础上推导得出的。

具体来说，Otsu算法首先基于直方图得到每个分割特性值的发生概率，并以阈值变量将分割特征值分为两类。

然后，算法计算出每一类的类内方差及类间方差，选取使类间方差最大或类内方差最小的阈值T作为最佳阈值。

当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好。

但总体来说，Otsu阈值最大类间方差法在图像分割、二值化等方面具有高效性，且该算法的稳定性能够得到保证。

图像的熵是一种特征的统计形式，它反映了图像中平均信息量的多少。

图像的一维熵表示图像中灰度分布的聚集特征所包含的信息量，令P i表示图像中灰度值为i的像素所占的比例，则定义灰度图象的一元灰度熵为：H= - E i=0255P i lnP i(其中，E i=0255表示从灰度0到255进行求和运算，因为公式无法输入，暂且这样表示)图象的一维熵可以表示图像灰度分布的聚集特征，却不能反映图像灰度分布的空间特征，为了表征这种空间特征，可以在一维熵的基础上引入能够反映灰度分布空间特征的特征量来组成图像的二维熵。

选择图像的邻域灰度均值作为灰度分布的空间特征量，与图像的像素灰度组成特征二元组，记为(i,j)，其中i表示像素的灰度值(0<=i<=255)，j表示领域灰度均值(0<=j<=255),P ij=(f(i,j))/N2即可反应某像素位置上的灰度值与其周围像素的灰度分布的综合特征，其中f(i, j)为特征二元组(i,j)出现的频数，N为图像的尺度，定义离散的图像二维熵为：H= - E i=0255E j=0255 (P ij lnP ij )依此构造的图像二维熵可以在反映图像所包含的信息量的前提下，突出反映图像中像素位置的灰度信息和像素邻域内灰度分布的综合特征；本文是Li和Lee关于一维最小交叉熵图像阈值分割的原文。

这里进行了简单翻译，有不足的地方请大家一起讨论完善。

后续有文章对信息熵学进行初窥，敬请期待。

摘要：通过最小化图像与其部分区域之间的交叉熵解决了图像分割中的阈值选取问题。

其中交叉熵基于两幅图像之间的像素运算得到，同时提出一种基于统计直方图的实现算法。

该方法提出了一种信息理论意义上针对二值图像的无偏估计，因而不需要了解图像灰度分布的先验知识。

1、介绍在众多的模式识别系统中，图像阈值化是重要的一个初步环节。

其中阈值的选取会影响到后期图像分割的精度以及效率。

本文的方法假设图像中的背景和待分割目标可以通过将图像的灰度值与选定的一个合适的阈值对比而区分。

OPENCV的二值化操作中，有一种“大津阈值处理”的方法，使用函数cvThreshold(image,image2,0,255,CV_THRESH_OTSU) 实现，该函数就会使用大律法OTSU得到的全局自适应阈值来进行二值化图片，而参数中的threshold不再起作用。

OTSU算法OTSU算法也称最大类间差法，有时也称之为大津算法，由大津于1979年提出，被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。

它是按图像的灰度特性,将图像分成背景和前景两部分。

因方差是灰度分布均匀性的一种度量,背景和前景之间的类间方差越大,说明构成图像的两部分的差别越大,当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。

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

设灰度图像灰度级是L，则灰度范围为[0,L-1]，利用OTSU算法计算图像的最佳阈值为：t = Max[w0(t) * (u0(t) - u)^2 + w1(t) * (u1(t) - u)^2)]其中的变量说明：当分割的阈值为t时，w0为背景比例，u0为背景均值，w1为前景比例，u1为前景均值，u为整幅图像的均值。

使以上表达式值最大的t，即为分割图像的最佳阈值。

以下是一段在OpenCV中实现的C语言程序，即一个使用OTSU算法提取图像阈值的函数，输入参数为一个图像指针，返回分割该图像的最佳阈值。

其中的变量说明：当分割的阈值为t时w0为背景像素点占整幅图像的比例u0为w0平均灰度w1为前景像素点占整幅图像的比例u1为w1平均灰度u为整幅图像的平均灰度公式:g = w0*pow((u-u0),2) + w1*pow((u-u1),2)int MyAutoFocusDll::otsuThreshold(IplImage *frame){const int GrayScale = 256;int width = frame->width;int height = frame->height;int pixelCount[GrayScale];float pixelPro[GrayScale];int i, j, pixelSum = width * height, threshold = 0;uchar* data = (uchar*)frame->imageData; //指向像素数据的指针for(i = 0; i < GrayScale; i++){pixelCount[i] = 0;pixelPro[i] = 0;}//统计灰度级中每个像素在整幅图像中的个数for(i = 0; i < height; i++){for(j = 0;j < width;j++){pixelCount[(int)data[i * width + j]]++; //将像素值作为计数数组的下标 }}//计算每个像素在整幅图像中的比例float maxPro = 0.0;int kk=0;for(i = 0; i < GrayScale; i++){pixelPro[i] = (float)pixelCount[i] / pixelSum;if( pixelPro[i] > maxPro ){maxPro = pixelPro[i];kk = i;}}//遍历灰度级[0,255]float w0, w1, u0tmp, u1tmp, u0, u1, u, deltaTmp, deltaMax = 0;for(i = 0; i < GrayScale; i++) // i作为阈值{w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;for(j = 0; j < GrayScale; j++){if(j <= i) //背景部分{w0 += pixelPro[j];u0tmp += j * pixelPro[j];}else //前景部分{w1 += pixelPro[j];u1tmp += j * pixelPro[j];}}u0 = u0tmp / w0;u1 = u1tmp / w1;u = u0tmp + u1tmp;deltaTmp = w0 * pow((u0 - u), 2) + w1 * pow((u1 - u), 2);if(deltaTmp > deltaMax){deltaMax = deltaTmp;threshold = i;}}return threshold;}接下来介绍OTSU方法的原理：********************************************************************* ********************************************************************* *********************OTSU法对于具有双峰性质的灰度图像或是彩色图像的某一通道的分割效果很好,程序为了增加健壮性加了个可以根据实际情况确定的修正值th_set.********************************************************************* ********************************************************************* ********************function y1=OTSU(image,th_set)image=imread('color1.bmp');gray=rgb2gray(image);%原图像的灰度图low_high=stretchlim(gray);%增强图像，似乎也不是一定需要gray=imadjust(gray,low_high,[]);% subplot(224);imshow(gray);title('after adjust');count=imhist(gray);[r,t]=size(gray);n=r*t;l=256;count=count/n;%各级灰度出现的概率for i=2:lif count(i)~=0st=i-1;breakendend%以上循环语句实现寻找出现概率不为0的最小灰度值for i=l:-1:1if count(i)~=0;nd=i-1;breakendend%实现找出出现概率不为0的最大灰度值f=count(st+1:nd+1);p=st;q=nd-st;%p和分别是灰度的起始和结束值u=0;for i=1:q;u=u+f(i)*(p+i-1);ua(i)=u;end%计算图像的平均灰度值for i=1:q;w(i)=sum(f(1:i));end%计算出选择不同k的时候，A区域的概率d=(u*w-ua).^2./(w.*(1-w));%求出不同k值时类间方差[y,tp]=max(d);%求出最大方差对应的灰度级th=tp+p;if th<th_setth=tp+p;elseth=th_set; %根据具体情况适当修正门限endy1=zeros(r,t);for i=1:rfor j=1:tx1(i,j)=double(gray(i,j));endendfor i=1:rfor j=1:tif (x1(i,j)>th)y1(i,j)=x1(i,j);elsey1(i,j)=0;endendend%上面一段代码实现分割% figure,imshow(y1);% title('灰度门限分割的图像');接下来介绍OTSU方法的原理以及用Ｃ语言实现：阈值将原图像分成前景、背景两个图像。