C++ OpenCV單峰三角閾值法Thresh_Unimodal詳解
需求說明
在對圖像進行處理時,經(jīng)常會有這類需求:想通過閾值對圖像進行二值化分割,以提取自己感興趣的區(qū)域,常見的閾值分割方法有常數(shù)分割、最大類間方差法、雙峰分割、三角法等等,不同的場景應用不同的閾值方法。
今天要講的方法,適合當圖像的直方圖具有明顯單峰特征時使用,結(jié)合了三角法的原理而設(shè)計,相比較OpenCV自帶的三角法,好處是可以根據(jù)自身需求合理修改函數(shù);如果用OpenCV庫的函數(shù),只有一個接口,若不能達到較理想的應用效果,就束手無策了。
下面介紹具體實現(xiàn)流程。
具體流程
1)取圖像的灰度圖,并遍歷統(tǒng)計0-255各個灰度值所出現(xiàn)的次數(shù)。
cv::Mat src = imread("test.jpg", 0);
cv::Mat hist = cv::Mat::zeros(1, 256, CV_32FC1);
for (int i = 0; i < src.rows; ++i)
{
for (int j = 0; j < src.cols; ++j)
{
hist.at<float>(0, src.at <uchar>(i, j))++;
}
}
2)去除0和255的直方圖數(shù)據(jù),這一步就是OpenCV三角法所沒有的。很多人可能不理解為什么要這一步,在你對圖像進行閾值化時如果提前進行了相關(guān)的運算,可能導致結(jié)果大于255的數(shù)值全部變?yōu)?55,或者數(shù)值低于0的數(shù)值全部變?yōu)?,這就使得0和255的數(shù)值其實涵蓋了許多數(shù)值,呈累加態(tài),很容易形成雙峰,這樣就很難找到我們真正想要的峰。例如0和255的數(shù)值都是10000左右,0略大一些,而我們的真峰是在250左右的灰度值,數(shù)值只有8000多,那么在后續(xù)閾值計算時就會因為峰的方向錯了而帶來毀滅性打擊。別覺得我說夸張了,只有自己去碰碰壁才能深刻領(lǐng)悟我說的。
hist.at<float>(0, 255) = 0; hist.at<float>(0, 0) = 0;
3)確認峰值位置,maxidx是峰值對應的灰度值,max是峰值高度,也是灰度值對應數(shù)據(jù)的個數(shù)。
float max = 0;
int maxidx = 0;
for (int i = 0; i < 256; ++i)
{
if (hist.at<float>(0, i) > max)
{
max = hist.at<float>(0, i);
maxidx = i;
}
}
4)判斷峰值在左側(cè)還是右側(cè),true為左側(cè),false為右側(cè)。
bool lr = maxidx < 127;
5)當在左側(cè)時,連接峰值(maxidx,max)和(255,0)點,用兩點建立直線公式,如下圖所示公式。 L的表達式可以轉(zhuǎn)換為Ax+By+C=0的形式,A是-max,B是maxidx-255,C是max*255,在結(jié)合距離公式可以計算出直方圖曲線上每個點到直線的距離,取距離最長的那個點作為閾值。

if (lr)
{
float A = float(-max);
float B = float(maxidx - 255);
float C = float(max * 255);
for (int i = maxidx + 1; i < 256; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
6)右側(cè)同理,連接峰值(maxidx,max)和(0,0)點,公式ABC如代碼所示。
else {
float A = float(-max);
float B = float(maxidx);
float C = 0.0f;
for (int i = 0; i < maxidx; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
7)二值化,完成。
result.setTo(255, src > maxdidx); idx = maxdidx; return result;
功能函數(shù)
// 單峰三角閾值法
cv::Mat Thresh_Unimodal(cv::Mat &src, int& idx)
{
cv::Mat result = cv::Mat::zeros(src.size(), CV_8UC1);
// 統(tǒng)計直方圖
cv::Mat hist = cv::Mat::zeros(1, 256, CV_32FC1);
for (int i = 0; i < src.rows; ++i)
{
for (int j = 0; j < src.cols; ++j)
{
hist.at<float>(0, src.at<uchar>(i, j))++;
}
}
hist.at<float>(0, 255) = 0;
hist.at<float>(0, 0) = 0;
// 搜索最大值位置
float max = 0;
int maxidx = 0;
for (int i = 0; i < 256; ++i)
{
if (hist.at<float>(0, i) > max)
{
max = hist.at<float>(0, i);
maxidx = i;
}
}
// 判斷最大點在哪一側(cè),true為左側(cè),false為右側(cè)
bool lr = maxidx < 127;
float maxd = 0;
int maxdidx = 0;
// 假設(shè)在左側(cè)
if (lr)
{
float A = float(-max);
float B = float(maxidx - 255);
float C = float(max * 255);
for (int i = maxidx + 1; i < 256; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
// 假設(shè)在右側(cè)
else {
float A = float(-max);
float B = float(maxidx);
float C = 0.0f;
for (int i = 0; i < maxidx; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
// 二值化
result.setTo(255, src > maxdidx);
idx = maxdidx;
return result;
}
C++測試代碼
#include <iostream>
#include <time.h>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
cv::Mat DrawHistImg(cv::Mat &hist);
cv::Mat Thresh_Unimodal(cv::Mat &src, int& idx);
int main()
{
cv::Mat src = imread("test.jpg", 0);
// 繪制均衡化后直方圖
cv::Mat hrI = DrawHistImg(src);
// 單峰三角閾值法
int thresh;
cv::Mat result = Thresh_Unimodal(src, thresh);
cout << " thresh: " << thresh << endl;
imshow("original", src);
imshow("hist", hrI);
imshow("result", result);
waitKey(0);
return 0;
}
// 繪制簡易直方圖
cv::Mat DrawHistImg(cv::Mat &src)
{
cv::Mat hist = cv::Mat::zeros(1, 256, CV_32FC1);
for (int i = 0; i < src.rows; ++i)
{
for (int j = 0; j < src.cols; ++j)
{
hist.at<float>(0, src.at <uchar>(i, j))++;
}
}
cv::Mat histImage = cv::Mat::zeros(540, 1020, CV_8UC1);
const int bins = 255;
double maxValue;
cv::Point2i maxLoc;
cv::minMaxLoc(hist, 0, &maxValue, 0, &maxLoc);
int scale = 4;
int histHeight = 540;
for (int i = 0; i < bins; i++)
{
float binValue = (hist.at<float>(i));
int height = cvRound(binValue * histHeight / maxValue);
cv::rectangle(histImage, cv::Point(i * scale, histHeight),
cv::Point((i + 1) * scale - 1, histHeight - height), cv::Scalar(255), -1);
}
return histImage;
}
// 單峰三角閾值法
cv::Mat Thresh_Unimodal(cv::Mat &src, int& idx)
{
cv::Mat result = cv::Mat::zeros(src.size(), CV_8UC1);
// 統(tǒng)計直方圖
cv::Mat hist = cv::Mat::zeros(1, 256, CV_32FC1);
for (int i = 0; i < src.rows; ++i)
{
for (int j = 0; j < src.cols; ++j)
{
hist.at<float>(0, src.at<uchar>(i, j))++;
}
}
hist.at<float>(0, 255) = 0;
hist.at<float>(0, 0) = 0;
// 搜索最大值位置
float max = 0;
int maxidx = 0;
for (int i = 0; i < 256; ++i)
{
if (hist.at<float>(0, i) > max)
{
max = hist.at<float>(0, i);
maxidx = i;
}
}
// 判斷最大點在哪一側(cè),true為左側(cè),false為右側(cè)
bool lr = maxidx < 127;
float maxd = 0;
int maxdidx = 0;
// 假設(shè)在左側(cè)
if (lr)
{
float A = float(-max);
float B = float(maxidx - 255);
float C = float(max * 255);
for (int i = maxidx + 1; i < 256; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
// 假設(shè)在右側(cè)
else {
float A = float(-max);
float B = float(maxidx);
float C = 0.0f;
for (int i = 0; i < maxidx; ++i)
{
float x0 = float(i);
float y0 = hist.at<float>(0, i);
float d = abs(A * x0 + B * y0 + C) / std::sqrt(A * A + B * B);
if (d > maxd)
{
maxd = d;
maxdidx = i;
}
}
}
// 二值化
result.setTo(255, src > maxdidx);
idx = maxdidx;
return result;
}
測試效果

圖1 原圖灰度圖

圖2 直方圖

圖3 閾值圖

圖4 閾值結(jié)果
通過imagewatch插件可以觀察閾值203是不是在距離最遠的位置,答案是肯定的。
如果函數(shù)有什么可以改進完善的地方,非常歡迎大家指出,一同進步何樂而不為呢~?
以上就是C++ OpenCV單峰三角閾值法Thresh_Unimodal詳解的詳細內(nèi)容,更多關(guān)于C++ OpenCV單峰三角閾值法的資料請關(guān)注腳本之家其它相關(guān)文章!
相關(guān)文章
C/C++判斷傳入的UTC時間是否當天的實現(xiàn)方法
在項目中經(jīng)常會顯示一個時間,如果這個時間在今日內(nèi)就顯示為時分秒,否則顯示為年月日,有需要的朋友可以參考一下2014-01-01

