具有掩膜功能的大津阈值法（Otsu）

学习记录…

前文说到大津阈值法是一种自适应的基于全局的阈值分割算法，只有在图像直方图分布为双峰的情况下才会呈现出一种比较好的分割效果，但是待分割图像直方图分布并不是每次都是理想的结果。可能会是光照的影响改变了原本为双峰的直方图分布，或者说背景本身就呈现出了两个灰度级，加上前景那就是三个灰度级了，等等一些情况都会造成使用Otsu分割失败。

试验用图如下：

如图所示：直接利用Otsu算法对原图进行分割，因为高灰度级分布的影响，我们想要的前景数字部分被归结为了背景，导致分割失败。期望分割图是利用上图中掩膜部分直方图而得到的。而在本实验中，直接用Otsu分割图取反形成mask掩膜（即对应于原图只有在掩膜图像中为白色的区域才参与运算）。

所以因为有掩膜的参与，我们可以实现任意形状区域的大津阈值法，不再是局限于全局图像。

其实方法本身不难，难点在你该通过什么方式才能得到合适的掩膜区域，而掩膜区域中必须包含待分割的前景。本实验中，采用先对图像进行阈值处理，从结果中提取掩膜，然后利用掩膜部分采用Otsu算法，得到正确分割图像。获取掩膜的方法多种，读者自行尝试。更多应用，欢迎一起交流学习！

试验代码如下，和前一节的内容基本上没有太大的出入：

int main()
{
	string path = "F:\\Speedlimit\\RGBCutImages10\\S825.jpg";

	Mat SrcImage = imread(path);
	if (!SrcImage.data) {
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}

	cv::Mat grayImage, OtsuImage, MaskImage;
	cvtColor(SrcImage, grayImage, COLOR_BGR2GRAY);
	threshold(grayImage, MaskImage, 0, 255, THRESH_BINARY_INV | THRESH_OTSU);
	
	//直方图参数定义
	cv::Mat gray_hist;
	const int histSize = 256;
	float range[] = { 0, 256 };
	const float* histRange[] = { range };
	bool uniform = true, accumulate = false;

	calcHist(&grayImage, 1, 0, MaskImage, gray_hist, 1, &histSize, histRange, uniform, accumulate);
	//获得归一化图像直方图
	cv::Mat norm_gray_hist = cv::Mat::zeros(gray_hist.size(), gray_hist.type());
	for (int i = 0; i < histSize; ++i)
		norm_gray_hist.at<float>(i) = gray_hist.at<float>(i) / cv::countNonZero(MaskImage);
	
	//全局均值和全局方差
	Mat mat_mean, mat_stddev;
	double gray_mean;

	meanStdDev(grayImage, mat_mean, mat_stddev, MaskImage);
	gray_mean = mat_mean.at<double>(0, 0);

	//遍历所有灰度，计算类间方差
	std::vector<double>sigma_ks(histSize);
	for (int k = 0; k < histSize; ++k)
	{
		double p1 = 0.0;				//p1类发生概率
		double m_k = 0.0;				//高达k阶累计平均灰度

		for (int i = 0; i <= k; ++i)
		{
			p1 += norm_gray_hist.at<float>(i);
			m_k += i * norm_gray_hist.at<float>(i);
		}
		if (p1 == 0.f || (1 - p1) == 0.f)					//分母不能为0
			sigma_ks[k] = 0.0;
		else
			sigma_ks[k] = (gray_mean * p1 - m_k) * (gray_mean * p1 - m_k) / (p1 * (1 - p1));
	}

	double max_Sigma_k = 0.0;
	std::vector<int>maxval_Ts;
	int Threshold_T = 0;									//最终输出的阈值T
	//找类间方差最大值
	for (int i = 0; i < sigma_ks.size(); ++i)
	{
		if (sigma_ks[i] > max_Sigma_k)
			max_Sigma_k = sigma_ks[i];
	}
	//找极大值对应的所有灰度值
	for (int i = 0; i < sigma_ks.size(); ++i)
	{
		if (max_Sigma_k == sigma_ks[i])
			maxval_Ts.push_back(i);
	}
	//如果极大值点不唯一，那么取对应各个极大值的各个k的平均值来得到最终阈值threshold_T
	for (int i = 0; i < maxval_Ts.size(); ++i)
		Threshold_T += maxval_Ts[i];
	Threshold_T = Threshold_T / maxval_Ts.size();
	//输出背景底色根据实际情况改为0或255
	cv::Mat dstImage(grayImage.size(), grayImage.type(), cv::Scalar(255));
	int width = grayImage.cols;
	int height = grayImage.rows;

	for (int i = 0; i < height; ++i)
	{
		for (int j = 0; j < width; ++j)
		{
			if (MaskImage.at<uchar>(i, j) == 255)
			{
				if (grayImage.at<uchar>(i, j) < Threshold_T)
					dstImage.at<uchar>(i, j) = 0;
			}
		}
	}

	imshow("src", SrcImage);
	cv::waitKey(0);
	return 0;
}

另外一张测试图片：

代码仅是一个小小功能的实现，测试的图像数量不是很多，有时间可以好好封装一下，后续不断更新一些经典的传统图像处理算法功能实现和改进，如有小伙伴需要类似的测试图片素材，可以私聊我。