【OpenCV3图像处理】阈值化操作得到二值图像（附Otsu算法和三角法原理与源码）

一、threshold()

函数原型：（opencv帮助文档）

     参数：

src	输入图像 (多通道或单通道, 8位或32位浮点型).
dst	输出图像（大小和类型和输入图像一致）
thresh	阈值
maxval	使用的thresh_binary和thresh_binary_inv阈值类型时的最大值
type	阈值化操作的类型 (see the cv::ThresholdTypes).

cv::ThresholdTypes可选的值及对应的函数方程如下：

函数说明：

函数使用固定的阈值分别应用于多通道。该函数通常用于从灰度图像中获取二值图像，或者用于去除噪声，即滤除过小或过大值的像素。函数支持的阈值有几种类型。它们由参数type决定。

当时用thresh_otsu或thresh_triangle，函数确定最优阈值使用Otsu（大津法）或三角形算法，而不是指定的阈值。函数返回Otsu或三角形算法计算得到的阈值。目前，使用大津法和三角法，输入必须为单通道8位的图像。

二、adaptiveThreshold( )

函数原型：

参数：

src	8位单通道图像
dst	相同大小的相同类型的目标图像
maxValue	分配给满足条件的像素的非零值，即满足添加的像素赋值为maxValue
adaptiveMethod	使用的自适应阈值分割算法, 可选值参考： cv::AdaptiveThresholdTypes
thresholdType	阈值类型，只有两个选项： THRESH_BINARY 或者THRESH_BINARY_INV，作用和threshold（）函数一致
blockSize	计算阈值时，每个像素采用的邻域的大小 ，填入的值：3, 5, 7,等等
C	常数减去平均值或加权平均数

cv::AdaptiveThresholdTypes：

ADAPTIVE_THRESH_MEAN_C	阈值T（x，y）是像素（x，y）的邻域的平均值减去常数C
ADAPTIVE_THRESH_GAUSSIAN_C	阈值T（x，y）是像素（x，y）的邻域的一个加权求和（与高斯窗互相关）减去常数C

三、示例

四、总结

我最喜欢使用的是threshold()配合参数thresh_otsu，也就是OTSU算法自动阈值化得到二值图。

附1：

    OTSU算法（大津法）是一种图像灰度自适应的阈值分割算法，是1979年由日本学者大津提出，并由他的名字命名的。大津法按照图像上灰度值的分布，将图像分成背景和前景两部分看待，前景就是我们要按照阈值分割出来的部分。背景和前景的分界值就是我们要求出的阈值。遍历不同的阈值，计算不同阈值下对应的背景和前景之间的类内方差，当类内方差取得极大值时，此时对应的阈值就是OTSU算法（大津法）所求的阈值。

Otsu实现思路

1. 计算0~255各灰阶对应的像素个数，保存至一个数组中，该数组下标是灰度值，保存内容是当前灰度值对应像素数

2. 计算背景图像的平均灰度、背景图像像素数所占比例

3. 计算前景图像的平均灰度、前景图像像素数所占比例

4. 遍历0~255各灰阶，计算并寻找类间方差极大值

C++代码实现，需要使用到opencv：

int OtsuAlgThreshold( Mat &image, Mat &binmask )
{
	if (image.channels() != 1)
	{
		cout << "Please input Gray-image!" << endl;
	}
	int T = 0; //Otsu算法阈值  
	double varValue = 0; //类间方差中间值保存
	double w0 = 0; //前景像素点数所占比例  
	double w1 = 0; //背景像素点数所占比例  
	double u0 = 0; //前景平均灰度  
	double u1 = 0; //背景平均灰度  
	double Histogram[256] = { 0 }; //灰度直方图，下标是灰度值，保存内容是灰度值对应的像素点总数  
 	uchar *data = image.data;

	double totalNum = 0.0; //像素总数
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			if ( binmask.at<uchar>(i, j) != 0 ) totalNum++;
		}
	}
	//计算灰度直方图分布，Histogram数组下标是灰度值，保存内容是灰度值对应像素点数  
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			if (binmask.at<uchar>(i, j) != 0) Histogram[data[i*image.step + j]]++;
		}
	}
	int minpos, maxpos;
	for (int i = 0; i < 255; i++)
	{
		if (Histogram[i] != 0)
		{
			minpos = i;
			break;
		}
	}
	for (int i = 255; i > 0; i--)
	{
		if (Histogram[i] != 0)
		{
			maxpos = i;
			break;
		}
	}

	for (int i = minpos; i <= maxpos; i++)
	{
		//每次遍历之前初始化各变量  
		w1 = 0;       u1 = 0;       w0 = 0;       u0 = 0;
		//***********背景各分量值计算**************************  
		for (int j = 0; j <= i; j++) //背景部分各值计算  
		{
			w1 += Histogram[j];   //背景部分像素点总数  
			u1 += j*Histogram[j]; //背景部分像素总灰度和  
		}
		if (w1 == 0) //背景部分像素点数为0时退出  
		{
			break;
		}
		u1 = u1 / w1; //背景像素平均灰度  
		w1 = w1 / totalNum; // 背景部分像素点数所占比例
							//***********背景各分量值计算**************************  

							//***********前景各分量值计算**************************  
		for (int k = i + 1; k < 255; k++)
		{
			w0 += Histogram[k];  //前景部分像素点总数  
			u0 += k*Histogram[k]; //前景部分像素总灰度和  
		}
		if (w0 == 0) //前景部分像素点数为0时退出  
		{
			break;
		}
		u0 = u0 / w0; //前景像素平均灰度  
		w0 = w0 / totalNum; // 前景部分像素点数所占比例  
							//***********前景各分量值计算**************************  

							//***********类间方差计算******************************  
		double varValueI = w0*w1*(u1 - u0)*(u1 - u0); //当前类间方差计算  
		if (varValue < varValueI)
		{
			varValue = varValueI;
			T = i;
		}
	}
	return T;
}

void OtsuAlgThreshold(Mat &image, Mat &binmask, Mat &imageOtsu)
{
	if (image.channels() != 1)
	{
		cout << "Please input Gray-image!" << endl;
	}
	int T = 0; //Otsu算法阈值  
	double varValue = 0; //类间方差中间值保存  
	double w0 = 0; //前景像素点数所占比例  
	double w1 = 0; //背景像素点数所占比例  
	double u0 = 0; //前景平均灰度  
	double u1 = 0; //背景平均灰度  
	double Histogram[256] = { 0 }; //灰度直方图，下标是灰度值，保存内容是灰度值对应的像素点总数  
	uchar *data = image.data;

	double totalNum = 0.0; //像素总数
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			int a = binmask.at<uchar>(i, j);
			if (a != 0) totalNum++;
		}
	}
	//计算灰度直方图分布，Histogram数组下标是灰度值，保存内容是灰度值对应像素点数  
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			if (binmask.at<uchar>(i, j) != 0) Histogram[data[i*image.step + j]]++;
		}
	}
	for (int i = 0; i < 255; i++)
	{
		//每次遍历之前初始化各变量  
		w1 = 0;       u1 = 0;       w0 = 0;       u0 = 0;
		//***********背景各分量值计算**************************  
		for (int j = 0; j <= i; j++) //背景部分各值计算  
		{
			w1 += Histogram[j];   //背景部分像素点总数  
			u1 += j*Histogram[j]; //背景部分像素总灰度和  
		}
		if (w1 == 0) //背景部分像素点数为0时退出  
		{
			break;
		}
		u1 = u1 / w1; //背景像素平均灰度  
		w1 = w1 / totalNum; // 背景部分像素点数所占比例  
		//***********背景各分量值计算**************************  

		//***********前景各分量值计算**************************  
		for (int k = i + 1; k < 255; k++)
		{
			w0 += Histogram[k];  //前景部分像素点总数  
			u0 += k*Histogram[k]; //前景部分像素总灰度和  
		}
		if (w0 == 0) //前景部分像素点数为0时退出  
		{
			break;
		}
		u0 = u0 / w0; //前景像素平均灰度  
		w0 = w0 / totalNum; // 前景部分像素点数所占比例  
		//***********前景各分量值计算**************************  

		//***********类间方差计算******************************  
		double varValueI = w0*w1*(u1 - u0)*(u1 - u0); //当前类间方差计算  
		if (varValue < varValueI)
		{
			varValue = varValueI;
			T = i;
		}
	}

	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			int a = binmask.at<uchar>(i, j);
			if (a != 0)
			{
				if (image.at<uchar>(i, j) > T)
				{
					imageOtsu.at<Vec3b>(i, j)[0] = 255;
					imageOtsu.at<Vec3b>(i, j)[1] = 0;
					imageOtsu.at<Vec3b>(i, j)[2] = 0;
				}
				else
				{
					imageOtsu.at<Vec3b>(i, j)[0] = 0;
					imageOtsu.at<Vec3b>(i, j)[1] = 0;
					imageOtsu.at<Vec3b>(i, j)[2] = 255;
				}
			}

		}
	}
}

附2：

三角法阈值分割：作者：gloomyfish