C#自动阈值选择：Otsu法编程实例(图像分割)

//Otsu阈值
        private void OtsuToolStripMenuItem_Click(object sender, EventArgs e)
        {
            if (m_Bitmap != null)
            {
                Rectangle rect = new Rectangle(0, 0, m_Bitmap.Width, m_Bitmap.Height);
                BitmapData bmpData = m_Bitmap.LockBits(rect, System.Drawing.Imaging.ImageLockMode.ReadWrite, m_Bitmap.PixelFormat);
                IntPtr ptr = bmpData.Scan0;
                int bytes = m_Bitmap.Width * m_Bitmap.Height;
                byte[] grayValues = new byte[bytes];
                System.Runtime.InteropServices.Marshal.Copy(ptr, grayValues, 0, bytes);

                byte temp = 0;               
                byte maxGray = 0;
                byte minGray = 255;
                int[] countPixel = new int[256];
                double mu1, mu2;

                //计算直方图
                for (int i = 0; i < grayValues.Length; i++)
                {
                    temp = grayValues[i];
                    countPixel[temp]++;
                    if (temp > maxGray)
                    {
                        //最大灰度等级
                        maxGray = temp;
                    }
                    if (temp < minGray)
                    {
                        //最小灰度等级
                        minGray = temp;
                    }
                }
                //Otsu法
                double w1 = 0, w2 = 0;
                double sum = 0;
                int numerator = 0;
                double sigma;
                double tempMax = 0;
                byte T = 0;
                for (int i = minGray; i < maxGray; i++)
                {
                    sum += i * countPixel[i];
                }
                for (int i = minGray; i < maxGray; i++)
                {
                    w1 += countPixel[i];
                    numerator += i * countPixel[i];
                    mu1 = numerator / w1;
                    w2 = grayValues.Length - w1;
                    mu2 = (sum - numerator) / w2;

                    //otsu公式
                    sigma = w1 * w2 * (mu1 - mu2) * (mu1 - mu2);

                    if (sigma>tempMax)
                    {
                        tempMax = sigma;
                        T = Convert.ToByte(i);
                    }
                }
                for (int i = 0; i < bytes; i++)
                {
                    if (grayValues[i] < T)
                    {
                        grayValues[i] = 0;
                    }
                    else
                    {
                        grayValues[i] = 255;
                    }
                }
                System.Runtime.InteropServices.Marshal.Copy(grayValues, 0, ptr, bytes);
                m_Bitmap.UnlockBits(bmpData);
                pictureBox1.Image = m_Bitmap;
            }
        }