最大类间方差法（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最大时，可以认为此时前景和背景差异最大，也就是此时的灰度是最佳阈值

 

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 = newint[256];               //图象直方图，共256个点

                 byte color;

                 byte* pline;

                 int n, n1, n2;

                 inttotal;                                           //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;

           }