基于OTSU(大津法)的图像分块的阈值分割

一、开发环境：

 Qt版本：Qt5.12.3
 VS版本：VS2017
 opencv版本：opencv-4.5.1-vc14_vc15

二、要求：实现基于图像分块+OTSU的图像分割

1.OTSU大津法实现

算法原理：OTSU的算法是假设存在阈值TH将所有图像分成两类C1，也叫前景（小于TH）和C2，也叫后景（大于TH），这两类的均值分别为m1和m2，整个图像的均值为mG。同时图像被分为C1类和C2类的概率分别为p1和p2。因此有：
$$m1\ast p1+m2\ast p2=mG;①$$
$$p1+p2=1;②$$
根据方差的知识，类间方差的表达式为:
$$\sigma ²=p1\ast (m1-mG)²+p2\ast (m2-mG)²;③$$
化简后，有
$$\sigma ²=p1\ast p2\ast (m2-m1)²;④$$
求能使σ²最大的灰度级最大的k为该图像的最佳阈值，这就是大津法(OTSU)。
实现：
输入：原图像（src）

// 大津法
double QtWidgetsApplication1::myOSTU(cv::Mat &src) {

	// 判断输入图像是否为空
	if (src.empty()) {
		//qDebug() << "empty";
		return -1;
	}
	// 图像转灰度图
	Mat dst;
	src.copyTo(dst);
	if (dst.channels() > 1) {
		cvtColor(dst, dst, COLOR_BGR2GRAY);
	}
	// 最佳阈值
	int mythreshold = 0;
	double maxVariance = 0;  // 最大类间方差 
	double p1 = 0, p2 = 0;  // 前景与背景像素点所占比例
	double m1 = 0, m2 = 0;  // 前景与背景像素值平均灰度
	double histogram[256] = { 0 };
	double sum = 0;  // 像素总和
	const double EPS = 1e-6;
	//统计256个bin，每个bin像素的个数
	int tmp = 0;
	int index = 0;
	for (int i = 0; i < dst.rows; i++) {
		for (int j = 0; j < dst.cols; j++) {
			index = i * dst.cols + j;
			tmp = (int)dst.data[index];
			histogram[tmp]++;
		}
	}
	// 计算所有像素点的频次总数
	for (int i = 0; i < 256; i++) {
		sum += histogram[i];
	}
	
	// 前景像素统计
	for (int i = 0; i <= 254; i++) {

		for (int j = 0; j <= i; j++) {
			p1 += histogram[j];//以i为阈值，统计前景像素个数
			m1 += j * histogram[j];//以i为阈值，统计前景像素灰度总和
		}
		p1 = p1 / sum; 
		m1 = m1 / p1;

		//背景像素统计
		for (int j = i + 1; j <= 255; j++) {
			m2 += j * histogram[j];//以i为阈值，统计背景像素灰度总和
			p2 += histogram[j];//以i为阈值，统计前景像素个数
		}
		p2 = 1 - p1;
		m2 = m2 / p2;
		double variance = p1 * p2*(m1 - m2)*(m1 - m2); //当前类间方差计算
		
		// 保留当前最佳阈值
		if (variance > maxVariance+EPS)
		{
			maxVariance = variance;
			mythreshold = i;
		
		}
		// 清零
		p1 = 0;
		p2 = 0;
		m1 = 0;
		m2 = 0;
	}

	// 返回最佳阈值
	return mythreshold;
}

2.图像分块+阈值分割

图像分块算法原理：在原图上需要的部分，然后分割
实现的方式有很多，本文的逻辑是
① 获取图像，对图像按比例进行切块，切块时应当注意不能整除的边界部分
② 用一个容器来保存存放这些小的图像
③ 获取需要处理的块数，为其最佳阈值申请堆空间
④ 调用上面写的大津法，得到每一小块的最佳阈值
⑤ 阈值分割，对小于最佳阈值的置为255（白），其他部分置为0（黑），这里看具体需要选择自己需要的灰度值
⑥ 图像合并
⑦ 显示图像
整个过程的源代码如下：
其中输入的参数分别是：
参数1：输入图像（src）
参数2：输出图像（dst）
参数3：阈值处理的算法(type – 大津法,这里是预留的接口，为后面添加其他的处理算法）
参数4：行需要切“几刀”（Rnum）
参数5：列需要切“几刀”（Cnum）

// 基于图像分块的阈值分割
void QtWidgetsApplication1::Threshold(cv::Mat &src, cv::Mat &dst, int type, int Rnum, int Cnum) {
	// 判断输入图像是否为空
	if (src.empty()) {
		return;
	}

	//图像分块
	QVector<Mat> ceilImg;
	
	Mat imageCut, roiImg; 
	Mat MergeImage(Size(src.cols, src.rows), src.type());  // 合并后图像
	int rowSize = 0;
	int colSize = 0;
	

	// 图像分块,每块大小为(src.cols / Cnum -- 列,src.rows / Rnum -- 行)
	for (int j = 0; j < Rnum; j++)
	{
		for (int i = 0; i < Cnum; i++)
		{
			Rect rect(i *(src.cols / Cnum), j * (src.rows / Rnum), src.cols / Cnum, src.rows / Rnum);

			imageCut = Mat(src, rect);
			roiImg = imageCut.clone();
			ceilImg.push_back(roiImg);
		}
	}
	
	// 行边界
	colSize = src.cols % Cnum;
	if (colSize > 0) {
		Rect rect(Cnum *(src.cols / Cnum), 0, colSize, src.rows);

		imageCut = Mat(src, rect);
		roiImg = imageCut.clone();
		ceilImg.push_back(roiImg);
	}
	// 列边界
	rowSize = src.rows%Rnum;
	if (rowSize > 0) {
		Rect rect(0, Rnum * (src.rows / Rnum), src.cols, rowSize);

		imageCut = Mat(src, rect);
		roiImg = imageCut.clone();
		ceilImg.push_back(roiImg);
	}
	// 大津法阈值分割
	 
	// 获取需要处理的图像块数
	int count = ceilImg.size();
	
	// 为每一块图像的最佳阈值申请堆空间
	double *myThreshold = new double[count];
	for (int t = 0; t < count; t++) {
		// 获取最佳阈值
		myThreshold[t] = myOSTU(ceilImg[t]);
		
		// 根据最佳阈值分割
		for (int i = 0; i < ceilImg[t].rows; i++) {
			for (int j = 0; j < ceilImg[t].cols; j++) {
				if (ceilImg[t].ptr<uchar>(i)[j] < myThreshold[t]) {
					ceilImg[t].ptr<uchar>(i)[j] = 0;
				}
				else {
					ceilImg[t].ptr<uchar>(i)[j] = 255;
				}
			}
		}
		//threshold(ceilImg[t], ceilImg[t], myThreshold[t], 255, THRESH_BINARY);
	}
	// 释放堆空间
	delete myThreshold;
	
	//图像合并
	int t = 0;

	for (int j = 0; j < Rnum; j++)
	{
		for (int i = 0; i < Cnum; i++)
		{
			Rect ROI(i *(src.cols / Cnum), j * (src.rows / Rnum), src.cols / Cnum, src.rows / Rnum);
			ceilImg[t].copyTo(MergeImage(ROI));
			t++;
		}
	}
	//合并边界
	if (colSize > 0) {
		Rect ROI(Cnum *(src.cols / Cnum), 0, colSize, src.rows);

		ceilImg[t].copyTo(MergeImage(ROI));
		t++;
	}
	if (rowSize > 0) {
		Rect ROI(0, Rnum * (src.rows / Rnum), src.cols, rowSize);
		ceilImg[t].copyTo(MergeImage(ROI));
		t++;
	}

	imshow("mergeimage", MergeImage);
	
}

3.效果图

3行1列的效果如下

3行2列的效果图如下

4行2列的效果如下

4行3列的效果图如下

由上面的结果图可以看出，这张图片分成4行2列，和3行2列的效果比较好
下面分享一下遇到的几个问题

1. 白线
   取矩形时应当注意，系统自带的优化问题，做运算时应当注意优先级。如图1：
   下面代码中i *(src.cols / Cnum) 写成 i *src.cols / Cnum，导致范围截断

Rect ROI(i *(src.cols / Cnum), j * (src.rows / Rnum), src.cols / Cnum, src.rows / Rnum);

2. 边缘问题 ---- 边界不做处理（一般边界重要信息不多）
   分为四种情况：行用R表示，列用C表示
   ① R、C均能整除,如下图
   

② R整除，C不能整除，如图

③ R不能整除，C整除，如图：

④ R、C均不能整除，如图：

初步解决，简单边界单独处理，如图：

3. imread(“xxx”,flag);问题
   当flag = 0,生成的灰度图与cvtColor(xx,xx,BGR转灰度)生成的灰度图有局部灰度不同。
   imread(“xxx”,0)效果图如图：

CvtColor效果图

• When using IMREAD_GRAYSCALE, the codec’s internal grayscale conversion will be used, if available. Results may differ to the output of cvtColor()
• 使用IMREAD_GRAYSCALE时，将使用编解码器的内部灰度转换（如果有）。结果可能与cvtColor（）的输出不同