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

一、开发环境：

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

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

1.OTSU大津法实现

算法原理：OTSU的算法是假设存在阈值TH将所有图像分成两类C1，也叫前景（小于TH）和C2，也叫后景（大于TH），这两类的均值分别为m1和m2，整个图像的均值为mG。同时图像被分为C1类和C2类的概率分别为p1和p2。因此有：
m 1 ∗ p 1 + m 2 ∗ p 2 = m G ; ① m1*p1+m2*p2 = mG;① m1∗p1+m2∗p2=mG;①
p 1 + p 2 = 1 ; ② p1+p2 = 1;② p1+p2=1;②
根据方差的知识，类间方差的表达式为:
σ ² = p 1 ∗ ( m 1 − m G ) ² + p 2 ∗ ( m 2 − m G ) ² ; ③ σ² = p1*(m1-mG)² + p2*(m2-mG)²;③ σ²=p1∗(m1−mG)²+p2∗(m2−mG)²;③
化简后，有
σ ² = p 1 ∗ p 2 ∗ ( m 2 − m 1 ) ² ; ④ σ² = p1* p2*(m2-m1)²;④ σ²=p1∗p2∗(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：
下面代码中i *(src.cols / Cnum) 写成 i *src.cols / Cnum，导致范围截断

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

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

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

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

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

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

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（）的输出不同