最大类间方差法(大津法OTSU)

        public static byte OstuThreshold(Bitmap bitmap)
        {
            byte threshold = 128;
            if (bitmap.PixelFormat == PixelFormat.Format8bppIndexed)
            {
                Rectangle rect = new Rectangle(0, 0, bitmap.Width, bitmap.Height);
                BitmapData bmpData = bitmap.LockBits(
                    rect,
                    ImageLockMode.ReadWrite,
                    bitmap.PixelFormat);

                int[] countPixel = new int[256];
                Array.Clear(countPixel, 0, 256);
                int offset = bmpData.Stride - bmpData.Width;
                int sum = bitmap.Height * bitmap.Width;
                double sumPixel = 0;
                double scale = 1 / (double)sum;

                unsafe
                {
                    byte* ptr = (byte*)bmpData.Scan0;
                    //calculate the count of each gray level
                    for (int y = 0; y < bmpData.Height; ++y)
                    {
                        for (int x = 0; x < bmpData.Width; ++x)
                        {
                            ++countPixel[*ptr];
                            ++ptr;
                        }
                        ptr += offset;
                    }
                }

                for (int i = 0; i < 256; ++i)
                {
                    sumPixel += i * countPixel[i];
                }

                int forePart = 0, backPart;
                double foreSum = 0, backSum;
                double foreAver, backAver;
                double minCost = -1, costFuc;

                for (int i = 0; i < 256; ++i)
                {
                    forePart += countPixel[i];
                    if (0 == forePart) continue;
                    backPart = sum - forePart;
                    if (0 == backPart) break;
                    foreSum += i * countPixel[i];
                    backSum = sumPixel - foreSum;
                    foreAver = foreSum / forePart;
                    backAver = backSum / backPart;
                    costFuc = forePart * scale * backPart * (foreAver - backAver) * (foreAver - backAver);
                    //find the cost function with the least value
                    if (costFuc > minCost)
                    {
                        minCost = costFuc;
                        threshold = (byte)i;
                    }
                }
                bitmap.UnlockBits(bmpData);
            }

            return threshold;
        }