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

 

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 = &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 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);
}




//
另外一个
//

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;
        }