多阈值处理利用双Otsu阈值-C#实现-基于EmguCv

对于一张图片，Otsu处理是寻找一个最适阈值进行分隔，而此算法是改进，寻找两个最适阈值，将图片分割成3个部分。

吃算法的理论可以参考《数字图像处理第三版-冈萨雷斯》10.3.6节，这里不予写出。

下面是代码：

   

     /// <summary>
        /// 最适双阈值处理
        /// </summary>
        /// <param name="image">要处理的图片</param>
        /// <param name="gray1">第一个区域要填充的灰度</param>
        /// <param name="gray2">第二个区域要填充的灰度</param>
        /// <param name="gray3">第三个区域要填充的灰度</param>
        /// <returns>返回要变化后的图片</returns>
        public static Mat DoubleThreshold(Mat image,int gray1=0,int gray2=127,int gray3=255)
        {
            Mat __mat = image.Clone();
            long[] _his = Histogram(__mat, 8);

            int _K1, _K2;//最适的两个阈值

            float _P1=0;//第一块类的的概率
            float _P2=0;//第二块类的的概率
            float _P3=0;//第三块类的的概率


            float _M1=0;//第一块类的的均值
            float _M2=0;//第二块类的的均值
            float _M3=0;//第三块类的的均值
            long _MN = _his.Sum();//总像素数
            float _MG;//图片平均灰度
          
            //用来存储最适阈值点和方差（因为point内部是存储两个int，所以直接拿来使用，不用自己再定义了）
            Dictionary<Point, double> _Variances = new Dictionary<Point, double>();

            //存储使方差最大的灰度集
            List<Point> _Ks= new List<Point>();

            #region 计算所有方差
            float _pi = 0;//表示一灰度像素所占的比例

            //分割点从零开始没什么意义所以从1开始
            for (int k=1;k<254;k++)
            {
                int k1;//表示一个分割的阈值
                int k2;//表示二个分割的阈值

                float _m1 = 0;//表示第一个区间的灰度累加均值
                float _m2 = 0;//表示第二个区间的灰度累加均值
                float _m3 = 0;//表示第三个区间的灰度累加均值

                _P1 = 0;
                _M1 = 0;
                _pi = 0;
                for(k1=0;k1<=k;k1++)
                {
                    _pi = (float)_his[k1] / _MN;
                    _P1 += _pi;
                    _m1 += k1 * _pi;
                }
                _M1 = _m1 / _P1;


               
             for(k2=k+1;k2<255;k2++)
                {

                    _pi = 0;
                    _M2 = 0;
                    _P2 = 0;
                    _m2 = 0;
                    for (int i=k+1;i<=k2;i++)
                    {
                        _pi = (float)_his[i] / _MN;
                        _P2 += _pi;
                        _m2 += i * _pi;
                    }
                    _M2 = _m2/ _P2;



                    _pi = 0;
                    _M3 = 0;
                    _P3 = 0;
                    _m3= 0;
                    for (int i = k2+1; i <= 255; i++)
                    {
                        _pi = (float)_his[i] / _MN;
                        _P3 += _pi;
                        _m3 += i * _pi;
                    }
                    _M3 = _m3 / _P3;


                    _MG = _m1 + _m2 + _m3;


                 var  _variance = _P1 * Math.Pow(_M1 - _MG, 2) + _P2 * Math.Pow(_M2 - _MG, 2) + _P3 * Math.Pow(_M3 - _MG, 2);

                    //这里将k1,j2减一是因为退出for循环时灰度k1,k2多加了一
                    _Variances.Add(new Point(k1-1, k2-1), _variance);
                }





            }
            #endregion

            double _MaxValue =   _Variances.Values.Max();
           foreach(var p in _Variances.Keys)
            {
                if (_Variances[p] == _MaxValue)
                {
                    _Ks.Add(p);
                }
            }

            _K1 = (int)_Ks.Average((x) => { return x.X; });
            _K2 = (int)_Ks.Average((x) => { return x.Y; });

            #region 像素操作
            unsafe
            {
                byte* data = (byte*)__mat.DataPointer.ToPointer();

                for (int row = 0; row < __mat.Height; row++)
                {
                    //data = data + row * image.Cols;
                    for (int col = 0; col < __mat.Width; col++)
                    {

                       if(0<=*data&&*data<_K1)
                        {
                            *data =(byte) gray1;
                        }
                        else if(_K1<=*data&&*data<_K2)
                        {
                            *data = (byte)gray2;
                        }
                        else
                        {
                            *data = (byte)gray3;
                        }

                        data++;
                    }
                }
            }
            #endregion
            return __mat;
        }

/// <summary>
        /// 图片直方图计算
        /// </summary>
        /// <param name="image">图片</param>
        /// <param name="depth">图片深度</param>
        /// <returns>放回直方图</returns>
        public static long[] Histogram(Mat image,int depth)
        {
            if (image.NumberOfChannels != 1)
            {
                throw new Exception("通道必须为1");
            }

            //提取直方图------------------------------------
            long[] _his = new long[256];
            for (int i = 0; i < 256; i++)
            {
                _his[i] = 0;
            }

            unsafe
            {
                byte* data = (byte*)image.DataPointer.ToPointer();

                for (int row = 0; row < image.Height; row++)
                {
                    //data = data + row * image.Cols;
                    for (int col = 0; col < image.Width; col++)
                    {
                        
                            _his[*data]++;
                       
                        data++;
                    }
                }
            }
            return _his;
        }

图片效果：                                   

原始图：                                                                                                                                                                     处理后：

原始图：                                                                                                                                                        处理后：