实用计算机视觉 -- 一种基于直方图的自动阈值计算方法

首先放上测试图片：

基于OpenCV的直方图显示代码：

Mat draw1DHistgoram(Mat &image, int channel){
float channel_range[] = { 0.0, 255.0 };
const float* channel_ranges = channel_range;
bool uniform = true; bool accumulate = false;
int histSize = 255;


// Draw the histograms
int hist_w = 512; int hist_h = 400;
int bin_w = cvRound((double)hist_w / histSize);
Mat display_image(hist_h, hist_w, CV_8UC3, Scalar(0));


if (channel == 1){
CV_Assert(image.channels() == 1);
Mat hist;
calcHist(&image, 1, 0, Mat(), hist,
1, &histSize, &channel_ranges, uniform, accumulate);
/// Normalize the result to [ 0, histImage.rows ]
normalize(hist, hist, 0, display_image.rows, NORM_MINMAX, -1, Mat());


for (int i = 1; i < histSize; ++i)
{
line(display_image, Point(bin_w*(i - 1), hist_h - cvRound(hist.at<float>(i - 1))),
Point(bin_w*(i), hist_h - cvRound(hist.at<float>(i))),
Scalar(0, 255, 255), 2, 8, 0);
}
}
else if (channel == 3){
CV_Assert(image.channels() == 3);
vector<Mat>bgr_plannes;
split(image, bgr_plannes);
Mat b_hist, g_hist, r_hist;
/// Compute the histograms:
calcHist(&bgr_plannes[0], 1, 0, Mat(), b_hist, 1, &histSize, &channel_ranges, uniform, accumulate);
calcHist(&bgr_plannes[1], 1, 0, Mat(), g_hist, 1, &histSize, &channel_ranges, uniform, accumulate);
calcHist(&bgr_plannes[2], 1, 0, Mat(), r_hist, 1, &histSize, &channel_ranges, uniform, accumulate);
/// Normalize the result to [ 0, histImage.rows ]
normalize(b_hist, b_hist, 0, display_image.rows, NORM_MINMAX, -1, Mat());
normalize(g_hist, g_hist, 0, display_image.rows, NORM_MINMAX, -1, Mat());
normalize(r_hist, r_hist, 0, display_image.rows, NORM_MINMAX, -1, Mat());
/// Draw for each channel
for (int i = 1; i < histSize; i++)
{
line(display_image, Point(bin_w*(i - 1), hist_h - cvRound(b_hist.at<float>(i - 1))),
Point(bin_w*(i), hist_h - cvRound(b_hist.at<float>(i))),
Scalar(255, 0, 0), 2, 8, 0);
line(display_image, Point(bin_w*(i - 1), hist_h - cvRound(g_hist.at<float>(i - 1))),
Point(bin_w*(i), hist_h - cvRound(g_hist.at<float>(i))),
Scalar(0, 255, 0), 2, 8, 0);
line(display_image, Point(bin_w*(i - 1), hist_h - cvRound(r_hist.at<float>(i - 1))),
Point(bin_w*(i), hist_h - cvRound(r_hist.at<float>(i))),
Scalar(0, 0, 255), 2, 8, 0);
}
}

return display_image;

}

灰度图片和RGB图片的直方图分别如下：

基于直方图的概率统计平均值是自动阈值计算的一种较好方法。

其中，p(g)是直方图概率分布，g是灰阶水平(0 ~ 255), Wb是背景像素比例，Wf是前景像素比例，当设置T的初始值为0且只进行一次迭代时，得到较为粗糙的自动阈值计算。

阈值化效果如下：

实现代码：

int computeOptimalTh(Mat &image, int init_threshold){
CV_Assert(image.channels() == 1 
&& init_threshold>=0 
&& init_threshold<=255);
float channel_range[] = { 0.0, 255.0 };
const float* channel_ranges = channel_range;
bool uniform = true; bool accumulate = false;
int histSize = 256;
Mat hist(histSize,1,CV_32FC1);
calcHist(&image, 1, 0, Mat(), hist,
1, &histSize, &channel_ranges, uniform, accumulate);
// Normalize the result to [ 0, histImage.rows ]
//normalize(hist, hist, 0, image.rows, NORM_MINMAX, -1, Mat());
Mat prob_hist(hist.size(), CV_32FC1);
cout << "HIST:" << hist.rows << " " << hist.cols << endl;

for (int row = 0; row < hist.rows; ++row){
prob_hist.at<float>(row, 0) = hist.at<float>(row, 0) / (image.rows*image.cols);
//cout << "probablity:" << prob_hist.at<float>(row, 0) << endl;
}
float Wb=0.0, Wf=0.0,Ub=0.0,Uf=0.0;
int T_pre, T_now;

T_pre = init_threshold;
for (int value = 0; value < init_threshold; ++value){
Wb += prob_hist.at<float>(value, 0);
}
Wf = 1.0 - Wb;
for (int value = 0; value < init_threshold; ++value){
Ub += (value * prob_hist.at<float>(value, 0)) / Wb;
//cout << "Ub:" << Ub << endl;
}
for (int value = init_threshold; value <= 255; ++value){
Uf += (value * prob_hist.at<float>(value, 0)) / Wf;
}


T_now = (int)(Ub + Uf) / 2;


return (int)T_now;
}

如果直方图近似为两个正态分布的叠加，便能在此算法的基础上进行全局最优阈值的自动计算。