全局阙值分割中的直方图算法和熵算法

前言：

阙值分割处理主要是根据灰度值信息提取前景，所以对前景物体与背景物体对比度较大的图像可以得到理想的分割效果，对于对比度较弱的图像，可以先提高图像的对比度，再进行阙值分割。

一、全局阙值分割

全局阙值分割是为所有的像素都设置同一个灰度值，大于或等于threshold的像素设为黑色或白色，相应地，小于threshold的像素设为白色或黑色。

在OpenCV3.2版本，重写了阙值函数threshold()，包含在"imgproc.hpp"头文件中，

double threshold( InputArray src, OutputArray dst,double thresh, double maxval, int type );enum ThresholdTypes {THRESH_BINARY   = 0, THRESH_BINARY_INV = 1, THRESH_TRUNC   = 2, THRESH_TOZERO   = 3, THRESH_TOZERO_INV = 4, THRESH_MASK = 7,THRESH_OTSU  = 8, THRESH_TRIANGLE  = 16
};

解释一下最后一个type，如果type为THRESH_BINARY，阙值分割规则为，当灰度大于阙值时，返回maxval，否则返回0；如果type为THRESH_BINARY_INV，则与前者规则相反。而如果type为THRESH_OTSU和THRESH_TRIANGLE，那么输入的图像矩阵只支持uchar类型，且thresh参数也作为输出参数，因为Ostu和TRIANGLE算法会自动计算thresh。

二、直方图技术进行阙值分割

当一幅图像有一个与背景呈现明显对比的物体的图像会具有包含双峰的直方图，而直方图技术则是根据这一“双峰”特点计算的。首先要找到这两个峰值，然后取两个峰值之间的波谷对应的灰度值，这个灰度值就是阙值，我们采用如下的算法来寻找直方图中的波谷。

假设输入图像为I，高位H，宽为W， $histogram_{I}$ 代表其对应的灰度直方图， $histogram_{I}(k)$ 代表灰度值等于k的像素点个数，其中 $0\leq k\leq 255$ ，

第一步：找到灰度直方图的第一个峰值，并找到其对应的灰度值也就是该峰值对应的位置，显然，第一个峰值就是灰度直方图的最大值，它对应的灰度值使用firstPeak表示。

第二步：找到灰度直方图的第二个峰值，并找到其对应的灰度值。注意：第二峰值不一定是直方图的第二大值，因为我们需要的峰值是希望在一定邻域内是最大值，如果第二大值和第一大值出现的很近，那么这种情况下第二大值就很明显不是一个峰值。我么使用如下公式计算最大峰值， $secondPeak=arg_{k}max\left \{ (k-firstPeak)^{2}*histogram_{I}(k) \right \}$ ，其中 $0\leq k\leq 255$ ，这个公式不仅考虑到了灰度值的大小，还将灰度值的位置与第一个峰值位置的距离作为参数。

第三步：找到了两个峰值后，就可以在两个峰值之间找到波谷灰度值，即该区域内灰度图的最小值的位置，然后设置为阙值进行阙值分割。其中，上面提到的阙值分割类型中的THRESH_TRIANGLE就和直方图技术类似。

实现代码：

//直方图技术自动计算全局阙值
int threshTwoPeaks(const Mat & image, Mat & threshImageOut) {Mat histogram = calcGrayHist(image);//计算第一个峰值位置Point firstPeakLoc;minMaxLoc(histogram,NULL,NULL,NULL,&firstPeakLoc);int firstPeak = firstPeakLoc.x;//计算第二个峰值位置Mat measureDist = Mat::zeros(Size(256,1),CV_32FC1);for (int k = 0; k < 256; k++) {int hist_k = histogram.at<int>(0,k);measureDist.at<float>(0, k) = pow((firstPeak - k), 2)*hist_k;}Point secondPeakLoc;minMaxLoc(measureDist,NULL,NULL,NULL,&secondPeakLoc);int secondPeak = secondPeakLoc.x;//得到波谷位置Point threshLoc;int thresh = 0;if (firstPeak < secondPeak) {minMaxLoc(histogram.colRange(firstPeak,secondPeak),NULL,NULL,&threshLoc);thresh = firstPeak + threshLoc.x+1;}else {minMaxLoc(histogram.colRange(secondPeak, firstPeak), NULL, NULL, &threshLoc);thresh = secondPeak + threshLoc.x + 1;}//根据波谷值进行阙值分割threshold(image,threshImageOut,thresh,255,THRESH_BINARY);return thresh;
}

注意：直方图技术只对灰度直方图有两个明显波峰的图像有较好的分割效果。

三、熵算法

信息熵(entropy)的概念来自信息论，假设信源信号u有n种取值，记为 $u_{1},u_{2}...u_{n}$ ，每一种信号出现的概率记为 $p_{1},p_{2}...p_{n}$ ，那么该信号的信号熵记为： $entroy(u)=-\sum_{i=1}^{n}p_{i}logp_{i}$ ，如果把图像看作是一个信源，假设输入图像为I， $normaHist_{I}$ 代表归一化的图像灰度直方图，那么对于8位图，就可以将其看为256个灰度符号，每一个符号出现的概率为 $normHist_{I}(k)$ ，

使用熵计算阙值的步骤如下：

第一步：计算归一化的灰度直方图I的累加概率直方图，又称零阶累积矩，记为 $cumuHist(k)=\sum_{i=0}^{k}normHist_{I}(i),k\in [0,255]$

所以，首先要获取灰度直方图，然后将其归一化得到normHist。

第二步：计算各个灰度级的熵，记为 $entropy(t)=-\sum_{k=0}^{t}normHist_{I}(k)\log (normHist_{I}(k))$ ，计算完之后将其存入一个1行256列的矩阵中

第三步：计算使 $f(t)=f_{1}(t)+f_{2}(t)$ 最大的t的值，即为阙值，其中

$f_{1}(t)=(entroy(t)/entroy(255)) * (log(normHist(t))/log\left \{ normHist(0),normHist(1),...,normHist(t) \right \})$

$f_{2}(t)=(1-entroy(t)/entroy(255)) * (log(1-normHist(t))/log\left \{ normHist(t+1),...,normHist(255) \right \})$

代码实现时，要特别注意在第二步和第三步的公式中，log函数要判断参数是否为0，进行相除操作时需要判断被除数是否为0，否则会出现异常。

//熵算法计算全局阙值
int threshByEntroy(const Mat & image, Mat & threshImageOut) {//获取灰度直方图Mat histogram = calcGrayHist(image);//归一化灰度直方图Mat normHist;histogram.convertTo(normHist,CV_32FC1,1.0/(image.rows*image.cols),0.0);//第一步，计算零阶累积矩阵Mat cumuHist=Mat::zeros(Size(256,1),CV_32FC1);for (int i = 0; i < 256; i++) {if (i == 0)cumuHist.at<float>(0, i) = normHist.at<float>(0, 0);elsecumuHist.at<float>(0, i) = cumuHist.at<float>(0,i-1)+normHist.at<float>(0,i);}//第二步，计算各个灰度级的熵Mat entroyHist = Mat::zeros(Size(256, 1), CV_32FC1);for (int i = 0; i < 256; i++) {float normHist_i=normHist.at<int>(0, i);if (i == 0) {if (normHist_i == 0)entroyHist.at<float>(0, i) = 0;elseentroyHist.at<float>(0, i) = -normHist_i*log10(normHist_i);}else {if (normHist_i == 0)entroyHist.at<float>(0, i) = entroyHist.at<float>(0, i - 1);elseentroyHist.at<float>(0, i) = entroyHist.at<float>(0, i - 1) - normHist_i*log10(normHist_i);}}printMat(entroyHist);//计算最大的f(t)，得到阙值Mat FHist = Mat::zeros(Size(256,1),CV_32FC1);float totalEntroy = entroyHist.at<float>(0, 255);for (int i = 0; i < 256; i++) {float cumuHist_i = cumuHist.at<float>(0, i);double maxVal1;minMaxLoc(normHist(Rect(0,0,i+1,1)),NULL,&maxVal1);float f1 = 0;if (cumuHist_i == 0 || cumuHist_i == 1 || maxVal1 == 0 || maxVal1 == 1 || totalEntroy == 0) {f1 = 0;}else {f1 = (entroyHist.at<float>(0, i) / totalEntroy)*(log10f(cumuHist_i) / log10f(maxVal1));}double maxVal2;minMaxLoc(normHist(Rect(i+1,0,255-i,1)),NULL,&maxVal2);float f2 = 0;if (cumuHist_i == 0 || cumuHist_i == 1 || maxVal2 == 0 || maxVal2 == 1) {f2 = 0;}else {if (totalEntroy == 0) {f2= log10f(1 - cumuHist_i) / log10f(maxVal2);}elsef2 = (1 - entroyHist.at<float>(0, i) / totalEntroy)*(log10f(1 - cumuHist_i) / log10f(maxVal2));}FHist.at<float>(0, i) = f1 + f2;}Point threshLoc;minMaxLoc(FHist,NULL,NULL,NULL,&threshLoc);//阙值分割操作threshold(image,threshImageOut,threshLoc.x,255,THRESH_BINARY);return threshLoc.x;
}