多阈值处理利用双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;

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

图片效果：                                   

原始图：                                                                                                                                                                     处理后：





原始图：                                                                                                                                                        处理后：