图像分割-基于Otsu算法的图像自适应阈值切割
1.Qtsu简介
Qtsu：Otsu算法（或称大津算法）。该算法由日本科学家大津展之（Nobuyuki Otsu）于1979年提出。这个算法看似简单，却与统计分析中的“方差分析”方法有非常深的渊源。我们今天在介绍该算法原理的基础之上，通过简单的C#代码来演示它的实现。我们在给灰度做二值化的时候，需要设定一个固定的阈值，我们并没有一个万能的阈值，而QSTU大津法则是根据灰度图像本身的信息，自动确定最佳阈值，实现以最佳阈值对灰度图进行二值化。需要注意的是，大津算法并不是直接进行二值化处理，而是通过统计学和概率知识计算计算出阈值，通过阈值再进行二值化。
一般情况下，我们对将一副图像分为背景和背景，前景一般是我们感兴趣的部分，背景一般是我们不感兴趣的部分。我们认为前景和背景之间的差异是比较大的，而前景部分中的像素是相似的，背景部分中的像素是相似的，同一类中差异小。不同类中差异大，如果存在一个阈值，使得图像分为了前景和背景，应当符合同一类中差异最小，不同类中差异最大，这个值就是最佳阈值。根据这个最佳阈值分割的图像效果应该是最好的（也就是将前景和背景彻底分割开）。
那么"差异"如何区分？其实就是用方差来表示：如果单个的数据其越偏离于中心，那么其方差值也就越大，如果同一类（前景或者背景）中方差越小，那么表示同一类中差异也小；不同类（前景和背景）之间方差越大，表示不同类中差异越大。
方差：其中μ表示期望，期望(Expectation, or expected value)是度量一个随机变量取值的集中位置或平均水平的最基本的数字特征；方差(Variance)是表示随机变量取值的分散性的一个数字特征。 方差越大，说明随机变量的取值分布越不均匀，变化性越强；方差越小，说明随机变量的取值越趋近于均值，即期望值。
对于离散型随机变量：
对于连续型随机变量：f(x)
2.Qtsu具体实现步骤：
（1）假设初始有个阈值T0，将图像分为两个部分，前景F和背景B。
（2）假设像素的总个数为N，前景像素个数为Nf，背景像素个数为Nb;
（3）假设图像的总灰度级为L-1,每个灰度级的像素个数为Ni，那么前景和背景像素总个数占总像素个数的概率分别满足如下的公式：
（4）前景和背景的灰度平均值分别为：
（5）那么整幅图像的灰度平均值为：
（6）那么，前景和背景之间的类方间差为：
（7）大津算法的目的就是求得一个阈值，使得第6步的类间方差最大，至于怎么求解，可以直接遍历或者其它优化算法。找出那个类间方差最大的阈值来分割即可。
下面是使用C#代码实现的Qtsu自适应阈值方法：
/// <summary> /// QTSU自适应阈值 /// </summary> /// <param name="mat">灰度图</param> /// <param name="minGray">统计的最小灰度</param> /// <param name="maxGray">统计的最大灰度</param> /// <returns></returns> int OtsuAlgThreshold(Mat mat, int minGray, int maxGray) { // 定义输出的目标直方图 Mat histogram = new Mat(); // 需要统计的通道索引 int[] channelIndex = new int[] { 0 }; // 每一维数值的统计范围 float[] ranges = new float[] { 0, 256 }; // 直方图横坐标的区间数。如果是10，则它会横坐标分为10份，然后统计每个区间的像素点总和。也就是bins，表示这个直方图分成多少份（即多少个直方柱） int[] hitSize = new int[] { 256 }; // 输入的Mat或Mat数组，并将灰度图赋值给它 VectorOfMat vMatImgs = new VectorOfMat(); vMatImgs.Push(mat); // 直方图计算 CvInvoke.CalcHist(vMatImgs, channelIndex, new Mat(), histogram, hitSize, ranges, false); // 找到直方图最大值 double minGrayValue = 0; double maxGrayValue = 0; System.Drawing.Point minLocationBlue = new System.Drawing.Point(); System.Drawing.Point maxLocationBlue = new System.Drawing.Point(); CvInvoke.MinMaxLoc(histogram, ref minGrayValue, ref maxGrayValue, ref minLocationBlue, ref maxLocationBlue); Image<Gray, int> histShowImage = histogram.ToImage<Gray, int>(); Dictionary<int, int> histdata = new Dictionary<int, int>();//灰度，灰度点个数 for (int i = 0; i < 256; i++) { histdata[i] = 0; } for (int i = minGray; i < maxGray; i++) { histdata[i] = histShowImage.Data[i, 0, 0]; } //去零大津法 int T = 0; //Otsu算法阈值 double varValue = 0; //类间方差中间值保存 double pb = 0;//背景像素点数占总点数的比例 double mb = 0;//背景所有像素点平均灰度 double pf = 0;//前景像素点数占总点数的比例 double mf = 0;//前景所有像素点平均灰度 double totalNum = mat.Rows * mat.Cols; //像素总数 for (int i = minGray; i <= maxGray; i++) { //每次遍历之前初始化各变量 pf = 0; mf = 0; pb = 0; mb = 0; //***********前景各分量值计算************************** for (int j = 0; j <= i; j++) //前景部分各值计算 { pf += histdata[j]; //前景部分像素点总数 mf += j * histdata[j]; //前景部分像素总灰度和 } if (pf == 0) //前景部分像素点数为0时退出 { continue; } mf = mf / pf; //前景像素平均灰度 pf = pf / totalNum; // 前景部分像素点数所占比例 //***********前景各分量值计算************************** //***********背景各分量值计算************************** for (int k = i + 1; k <=maxGray; k++) { pb += histdata[k]; //背景部分像素点总数 mb += k * histdata[k]; //背景部分像素总灰度和 } if (pb == 0) //背景部分像素点数为0时退出 { break; } mb = mb / pb; //背景像素平均灰度 pb = pb / totalNum; // 背景部分像素点数所占比例 //***********背景各分量值计算************************** //***********类间方差计算****************************** double varValueI = pb * pf * (mf - mb) * (mf - mb); //当前类间方差计算 if (varValue < varValueI) { varValue = varValueI; T = i; } } return T; }
