带有 mask 的 OTSU 自适应阈值

上一篇 图像二值化方法--OTSU（最大类间方差法） 留下一个问题：使用库函数 

threshold(src, dst, thresh, maxVal, CV_THRESH_OTSU)

来计算阈值，冒出了一个不合群的阈值，原因未知。

分析灰度图，发现光照亮度不均，中间较亮，四周较暗（进行直方图均衡化后更直观）。想把四周较暗的区域作为模板，计算该区域的OTSU阈值，作为整个图像的分割阈值。即把较暗区域的像素点和黑色目标区域的像素点作为对象求其OTSU阈值，这样计算出来的阈值更低。因此，就需要一个带 mask 的计算OTSU阈值的函数。这个想法是学习了以下内容后产生的，该作者改写了原函数。（第一次看到改写原函数，十分感谢原文作者）

opencv threshold with mask

thresh = getThreshVal_Otsu_8u(src); 该函数源码在  opencv\sources\modules\imgproc\src\thresh.cpp 中，文章开头的链接中也贴出了。

改写后的函数：

double getThreshVal_Otsu_8u_mask(const Mat src, const Mat& mask){const int N = 256;int M = 0;int i, j, h[N]={0};for( i=0; i<src.rows ; i++ ){const uchar* psrc = src.ptr(i);const uchar* pmask = mask.ptr(i);for( j=0; j<src.cols ; j++ ){if( pmask[j] ){h[ psrc[j] ]++;++M;}}}double mu =0, scale = 1. /(M);for(i=0; i<N; i++)mu += i * (double)h[i];mu *= scale;double mu1 = 0, q1 = 0;double max_sigma = 0, max_val = 0;for(i=0; i<N; i++){double p_i, q2, mu2, sigma;p_i = h[i] * scale;mu1 *= q1;q1 += p_i;q2 = 1. -q1;if(std::min(q1, q2) < FLT_EPSILON || std::max(q1, q2) > 1. - FLT_EPSILON)continue;mu1 = (mu1 + i*p_i) / q1;mu2 = (mu - q1*mu1) / q2;sigma = q1 * q2 * (mu1-mu2)*(mu1-mu2);if( sigma > max_sigma ){max_sigma = sigma;max_val = i;}}return max_val;}

main.cpp 

主要步骤见注释

int main(){// Load an imageMat src = imread("a (16).bmp");if(src.empty())cout<<"src is empty"<<endl;Mat imgLab;cvtColor(src, imgLab, CV_BGR2Lab);vector<Mat> imgLabChannels;split(imgLab, imgLabChannels);Mat imgLab_L = imgLabChannels[0];//Mat imgLab_a = imgLabChannels[1];//Mat imgLab_b = imgLabChannels[2];imshow("imgLab_L", imgLab_L); imwrite("imgLab_L.bmp", imgLab_L);// 直接使用OTSUMat imgThre;double thresh = threshold(imgLab_L, imgThre, 1, 255, CV_THRESH_OTSU);cout << "thresh=" << thresh << endl ;imshow("imgThre", imgThre); imwrite("imgThre.bmp", imgThre);//imgLab_b直方图均衡化Mat  imgEhi;equalizeHist(imgLab_L, imgEhi);imshow("imgEhi", imgEhi); imwrite("imgEhi.bmp", imgEhi);// 直接使用OTSUMat imgEhiThre;double thresh_Ehi = threshold(imgEhi, imgEhiThre, 1, 255, CV_THRESH_OTSU);imshow("imgEhiThre", imgEhiThre); imwrite("imgEhiThre.bmp", imgEhiThre);// 提取出其中光线较亮的中间部分: mask2// closeMat mask1(src.size(), CV_8UC1, Scalar(0));Mat elememt1 = getStructuringElement(MORPH_RECT, Size(3,3));morphologyEx(imgEhiThre, mask1, MORPH_DILATE, elememt1, Point(-1,-1), 3);//morphologyEx(mask1, mask1, MORPH_ERODE, elememt1, Point(-1,-1), 2);//边界置0，排除干扰mask1.rowRange(0,5).setTo(0);mask1.rowRange(src.rows-2, src.rows).setTo(0);mask1.colRange(0,2).setTo(0);mask1.colRange(src.cols-2, src.cols).setTo(0);imshow("mask1", mask1); imwrite("mask1.bmp", mask1);// 筛选连通域，把中间较亮区域提取出来，并且把该区域内的孔洞填充vector<vector<Point>> contours_mask1;vector<Vec4i> hierarchy_mask1;findContours(mask1, contours_mask1, hierarchy_mask1, CV_RETR_LIST, CV_CHAIN_APPROX_SIMPLE);Mat mask2(src.size(), CV_8UC1, Scalar(0));for(int i=0; i<contours_mask1.size(); ++i){double area = contourArea(contours_mask1[i]);if(area>100000)drawContours(mask2, contours_mask1, i, 255, CV_FILLED, 8, hierarchy_mask1);}imshow("mask2", mask2); imwrite("mask2.bmp", mask2);    Mat mask2_inv = 255. - mask2;imshow("mask2_inv", mask2_inv); imwrite("mask2_inv.bmp", mask2_inv);    // 作为模板// 使用带 mask 的OTSU阈值Mat mask3(src.size(), CV_8UC1, Scalar(0));double thresh_b = getThreshVal_Otsu_8u_mask(imgLab_L, mask2_inv);threshold(imgLab_L, mask3, thresh_b+40, 255, CV_THRESH_BINARY_INV);//threshold_mask(imgEhi, mask3, thresh_w, 255, CV_THRESH_OTSU, mask2);cout << "thresh_b=" << thresh_b << endl ;imshow("mask3", mask3); imwrite("mask3.bmp", mask3);    // 作为模板waitKey(0);return 0;}

（参考文章开头链接处的文章，内容是延续的，测试图也是该文的）

一些中间结果图：分别为 imgEhiThre, mask2

  

阈值对比结果：

从结果中看，这种方法算得的阈值基本统一。

但是只测试了这四张图，测试更多的图可能会发现更多问题。