基於sobel算子的邊緣檢測算法的C++實現及改進
實驗內容
- 對已有的sobel算子邊緣檢測算法的學習
- 對已有算法進行C++實現
- 對已有算法進行改進
已有的邊緣檢測算法流程
- 對圖像進行灰度化得到灰度圖像
- 對灰度圖像進行高斯模糊去除部分噪聲
- 利用sobel算子作為卷積核對圖像進行卷積，計算梯度
- 用非極大值抑制法進行邊緣細化
- 利用雙閾值進行噪聲去除
算法實現
圖像灰度化
cv::Mat Graying(cv::Mat img) {
int row = img.rows;
int col = img.cols;
// cv::imshow("灰度圖", img);
cv::Mat grayImg = Mat(row, col, CV_8U);
cout << row << '\n' << col;
for(register int i = 0; i < row; ++i) {
Vec3b *p = img.ptr<Vec3b>(i);
uchar *p2 = grayImg.ptr<uchar>(i);
for(register int j = 0; j <=col; ++j) {
*p2 = (*p)[0] * 0.114 + (*p)[1] * 0.587 + (*p)[2] * 0.299;
//Gray = 0.2989*R+0.5870*G+0.1140*B 注意opencv的Mat類型是BGR順序
p2++;
p++;
}
}
return grayImg;
}
灰度化實現效果如下:
原始圖像
灰度圖像
邊緣檢測主體實現
高斯模糊
GaussianBlur(grayImg, grayImg, Size(5, 5), BORDER_DEFAULT); //進行5*5的高斯模糊（這個是最后懶了，覺得代碼只是繁瑣就不想手寫，調的庫
sobel算子濾波
sobel算子分為x方向和y方向，在某個像素點處分別用兩個算子卷積求出x方向和y方向的導數值，再通過加權平均可以求出在某個像素點處的梯度值，而由於圖像邊緣就是圖像灰度值發生劇烈變化的地方，梯度值大的地方是邊緣的可能就大，因此sobel算子可以檢測出邊緣
\(G_x = \left[ \begin{matrix} -1 & 0 & 1\\-2 & 0 & 2\\-1 & 0 & 1 \end{matrix} \right] * A\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,G_y = \left[ \begin{matrix} -1 & -2 & -1\\0 & 0 & 0\\1 & 2 & 1 \end{matrix} \right] * A\\\\G = \sqrt{{G_x}^2+{G_y}^2}\)
cv::Mat filterlengthways(cv::Mat img1) {
cv::Mat img;
img1.copyTo(img);
int row = img.rows;
int col = img.cols;
vector<int> temp[row];
for (int i = 0; i < row; ++i) {
uchar *p = img.ptr<uchar>(i);
for (int j = 0; j < col; ++j) {
temp[i].push_back((int) (*p));
p++;
}
}
cv::Mat zero = Mat(cv::Size(col, row), CV_8UC1, Scalar(0));
zero.copyTo(img);
for (register int i = 1; i < row - 1; ++i) { //豎直方向濾波，算子為
for (register int j = 1; j < col - 1; ++j) { //-1 0 1
int temp1 = 0; //-2 0 2
temp1 -= temp[i - 1][j - 1]; //-1 0 1 該算子可以求出水平方向的導數
temp1 -= temp[i][j - 1] * 2;
temp1 -= temp[i + 1][j - 1];
temp1 += temp[i - 1][j + 1];
temp1 += temp[i][j + 1] * 2;
temp1 += temp[i + 1][j + 1];
if (temp1 < 0) temp1 = -temp1;
if (temp1 > 255) temp1 = 255;
img.at<uchar>(i, j) = temp1;
}
}
return img;
}
cv::Mat filtercrossways(cv::Mat img1) {
cv::Mat img;
img1.copyTo(img);
int row = img.rows;
int col = img.cols;
vector<int> temp[row];
for (int i = 0; i < row; ++i) {
uchar *p = img.ptr<uchar>(i);
for (int j = 0; j < col; ++j) {
temp[i].push_back((int) (*p));
p++;
}
}
cv::Mat zero = Mat(cv::Size(col, row), CV_8UC1, Scalar(0));
zero.copyTo(img);
for (register int i = 1; i < row - 1; ++i) {
for (register int j = 1; j < col - 1; ++j) { //水平方向濾波, 算子為
int temp1 = 0; //-1 -2 -1
temp1 -= temp[i - 1][j - 1]; // 0 0 0
temp1 -= temp[i - 1][j] * 2; // 1 2 1 求出豎直方向的導數
temp1 -= temp[i - 1][j + 1];
temp1 += temp[i + 1][j - 1];
temp1 += temp[i + 1][j] * 2;
temp1 += temp[i + 1][j + 1];
if (temp1 < 0) temp1 = -temp1;
if (temp1 > 255) temp1 = 255;
img.at<uchar>(i, j) = temp1;
}
// puts("");
}
return img;
}
vector<double> gradient;//求出梯度方向和大小，並進行梯度方向的規范化
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
int a = img1.at<uchar>(i, j), b = img2.at<uchar>(i, j);
int temp = sqrt(a * a + b * b);//幾何平均求出梯度大小
gradient.push_back(atan2(b, a));
if (gradient[getNum(i, j, col)] < 0) gradient[getNum(i, j, col)] += 2 * pi;//使得梯度方向為0到2pi，方便處理
img.at<uchar>(i, j) = temp;
}
}
豎直方向濾波
水平方向濾波
幾何平均求出梯度大小
非極大值抑制
由於sobel算子計算出的邊緣過於模糊和粗大，因此需要進行邊緣的細化。
梯度
不妨假設圖像邊緣部分的梯度如圖，且之前的sobel算子濾波會將橙線之上的部分都當成邊緣，因此邊緣會粗大而模糊
於是就產生一種想法即只保留頂點所對應的像素點作為邊緣，其它部分置為0，這就產生了非極大值抑制
之前sobel算子濾波時計算出了每個像素點梯度的大小和方向
\[G = \sqrt{{G_x}^2+{G_y}^2} \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\arctan({G_y \over G_x}) \]
因為圖像中方格是離散的，因此對於梯度的方向進行區域的划分，如圖所示
區域划分
對於每個像素點，計算出梯度大小和梯度方向后，將其梯度大小與對應區域相鄰像素的梯度大小進行比較，如果該像素點的梯度大小比任意一個對應區域相鄰像素小，則該像素點置為0，反之保持原樣，這樣就實現了邊緣的細化
const double pi = acos(-1);
const double zone[8] = {atan2(1, 3), atan2(3, 1), atan2(3, -1), atan2(1, -3), 2 * pi + atan2(-1, -3), 2 * pi + atan2(-3, -1), 2 * pi + atan2(-3, 1), 2 * pi + atan2(-1, 3)};
//進行區域划分，實際上是4個，但為了編碼方便，化成8個
//注意之前求梯度方向時對atan2的結果都加上了2*pi
int getType(double t) { //返回梯度方向所屬區域
for (int i = 0; i < 8; ++i) {
if (t >= zone[7] || t <= zone[0]) return 0;
if (t >= zone[0] && t <= zone[1]) return 1;
if (t >= zone[1] && t <= zone[2]) return 2;
if (t >= zone[2] && t <= zone[3]) return 3;
if (t >= zone[3] && t <= zone[4]) return 0;
if (t >= zone[4] && t <= zone[5]) return 1;
if (t >= zone[5] && t <= zone[6]) return 2;
if (t >= zone[6] && t <= zone[7]) return 3;
}
throw "error";
}
void nonMaximumSuppression(Mat img, vector<double> &gradient) {
int col = img.cols, row = img.rows;
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
switch (getType(gradient[getNum(i, j, col)])) {
case 0:
if (check(i, j - 1, row, col) &&
img.at<uchar>(i, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i, j + 1, row, col) &&
img.at<uchar>(i, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 1:
if (check(i - 1, j + 1, row, col) &&
img.at<uchar>(i - 1, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i + 1, j - 1, row, col) &&
img.at<uchar>(i + 1, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 2:
if (check(i - 1, j, row, col) &&
img.at<uchar>(i - 1, j) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i + 1, j, row, col) &&
img.at<uchar>(i + 1, j) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
case 3:
if (check(i + 1, j + 1, row, col) &&
img.at<uchar>(i + 1, j + 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
if (check(i - 1, j - 1, row, col) &&
img.at<uchar>(i - 1, j - 1) > img.at<uchar>(i, j)) {
img.at<uchar>(i, j) = 0;
}
break;
}
}
}
}
非極大值抑制后的邊緣
雙閾值判定
非極大值抑制后的圖像仍存在許多噪聲，可以將灰度值小於某個固定閾值的像素點置為0，可以一定程度上消除噪聲
void imageBinaryzation(Mat img, int thresholdVal = 0) { //嘗試對邊緣進行單固定閾值，但效果不好
Mat img1;
img.copyTo(img1);
Mat_<uchar>::iterator it1 = img1.begin<uchar>();
Mat_<uchar>::iterator itend1 = img1.end<uchar>();
for (; it1 != itend1; ++it1) {
(*it1) > thresholdVal ? 0 : (*it1) = 0; //如果像素點灰度值大於閾值，灰度不變，否則灰度置為0;
}
OutPut("Binaryzation", img1);
}
固定單閾值
然而固定單閾值法有許多問題
- 對於不同類型的圖像，不同要求的邊緣檢測，閾值不同，得到效果好的閾值很麻煩
- 極容易造成邊緣的中斷
因此采取雙閾值法可以一定程度上解決問題
我們定義高閾值和低閾值，對於灰度值高於高閾值的像素點我們將其確定為邊緣，而灰度值低於低閾值的像素點我們認為它一定不是邊緣，對於灰度值介於高閾值和低閾值之間的像素點，它有可能是邊緣，當它相鄰8個像素點存在確定的邊緣時，它也會變成確定的邊緣，否則它被判定為不是邊緣它的灰度值將被置為0
void doubleThresholdDetection(Mat &img1) {
int col = img1.cols, row = img1.rows;
cv::Mat temp;
img1.copyTo(temp);
for (int i = 0; i < row; ++i) {
for (int j = 0; j < col; ++j) {
if (temp.at<uchar>(i, j) > highTh) {
continue;
} else {
if (temp.at<uchar>(i, j) < lowTh) {
img1.at<uchar>(i, j) = 0;
continue;
} else {
if (lowthcheck(temp, i - 1, j - 1, row, col) ||
lowthcheck(temp, i - 1, j, row, col) ||
lowthcheck(temp, i - 1, j + 1, row, col) ||
lowthcheck(temp, i, j - 1, row, col) ||
// lowthcheck(temp, i, j , row, col) ||
lowthcheck(temp, i, j + 1, row, col) ||
lowthcheck(temp, i + 1, j - 1, row, col) ||
lowthcheck(temp, i + 1, j, row, col) ||
lowthcheck(temp, i + 1, j + 1, row, col)
) {
temp.at<uchar>(i, j) = highTh + 1;
continue;
} else {
img1.at<uchar>(i, j) = 0;
}
}
}
}
}
}
雙閾值法
可以很明顯的發現，雙閾值法較好地解決了單閾值法會造成邊緣中斷的問題，並且由於介於高低閾值之間的點仍能被保留，因此檢測效果較好的閾值范圍變大，閾值更容易調整
改進的邊緣細化方法
由於非極大值抑制對邊緣進行細化時，只會保留梯度最大的一個點，因此邊緣大多被細化成了一個像素點寬度的線，並且有可能使較粗的邊緣變成兩條線，這就使得邊緣鋸齒化嚴重結果不夠美觀，因此我試圖用另一種思路進行邊緣細化
改進的細化效果
非極大值抑制細化效果
\[sumLengthWay = \left[ \begin{matrix} 0 & 0 & 0 &0 & 2 & 0 &0 & 0 & 0\\ 0 & 0 & 0 &2 & 2 & 2 &0 & 0 & 0\\ 0 & 0 & 0 &2 & 2 & 2 &0 & 0 & 0\\ 0 & 0 & 1 &1 & 1 & 1 &1 & 0 & 0\\ 0 & 0 & 1 &1 & 1 & 1 &1 & 0 & 0\\ 0 & 0 & 1 &1 & 1 & 1 &1 & 0 & 0\\ 0 & 0 & 0 &2 & 2 & 2 &0 & 0 & 0\\ 0 & 0 & 0 &2 & 2 & 2 &0 & 0 & 0\\ 0 & 0 & 0 &0 & 2 & 0 &0 & 0 & 0\\ \end{matrix} \right] * A\]
\[sumCrossWay = \left[ \begin{matrix} 0 & 0 & 0 &0 & 0 & 0 &0 & 0 & 0\\ 0 & 0 & 0 &0 & 0 & 0 &0 & 0 & 0\\ 0 & 0 & 0 &1 & 1 & 1 &0 & 0 & 0\\ 0 & 2 & 1 &1 & 1 & 1 &1 & 2 & 0\\ 2 & 2 & 1 &1 & 1 & 1 &1 & 2 & 2\\ 0 & 2 & 1 &1 & 1 & 1 &1 & 2 & 0\\ 0 & 0 & 0 &1 & 1 & 1 &0 & 0 & 0\\ 0 & 0 & 0 &0 & 0 & 0 &0 & 0 & 0\\ 0 & 0 & 0 &0 & 0 & 0 &0 & 0 & 0\\ \end{matrix} \right] * A\]
\(cntCrossWay\) 和 \(cntLengthWay\) 分別為上述兩個矩陣的系數和
\(averCrossWay = {sumCrossWay \over cntCrossWay}\)
\(averLengthWay = {sumLengthWay \over cntLengthWay}\)
\(threshold Value = 0.92 *\sqrt{averLengthWay ^2+averCrossWay^2}\)
代碼實現:
void autoImageBinaryzation6(Mat img) {
int row = img.rows, col = img.cols;
Mat img1;
img.copyTo(img1);
int averCrossWay, averLengthWay, sumCrossWay, sumLengthWay, cntCrossWay, cntLengthWay, aver;
for(int j = 0; j < col; ++j) {
for(int i = 0; i < row; ++i) {
sumCrossWay = sumLengthWay = cntCrossWay = cntLengthWay = 0;
for(int k = -4; k <= 4; ++k) {
for(int g = (k <= 1 ? -2 : ((k <= 3 ) ? -1 : 0)) ; g <= (k <= 1 ? 2 : ((k <=3) ? 1 : 0)); ++g) {
if(check(i + g, j + k, row, col)) {
sumLengthWay += img.at<uchar>(i + g, j + k) * (1 + (abs(k) >= 2));
cntLengthWay += (1 + (abs(k) >= 2));
}
if(check(i + k, j + g, row, col)) {
sumCrossWay += img.at<uchar>(i + k, j + g) * (1 + (abs(k) >= 2));
cntCrossWay += (1 + (abs(k) >= 2));
}
}
}
averLengthWay = sumLengthWay / cntLengthWay;
averCrossWay = sumCrossWay / cntCrossWay;
aver = sqrt(averCrossWay * averCrossWay + averLengthWay * averLengthWay);
img1.at<uchar>(i, j) > aver * 0.92 ? 0 : img1.at<uchar>(i, j) = 0;
}
}
OutPut("autoBinaryzation6", img1);
doubleThresholdDetection(img1);
OutPut("autoBinaryzation6dtd", img1);
OutPut("myfilter", img1);
}
該方法看起來時間復雜度略有增大，但averLengthWay和averCrossWay可以進行動態維護，不需要每次進行遍歷求值，代碼實現中直接for循環遍歷求值只是為了呈現更加直觀。 且只涉及整型操作，因此實際運行速度並不比涉及到浮點運算以及atan2操作的非極大值抑制有劣勢。
雙閾值檢測后
效果展示
細化前
改進的細化效果
非極大值抑制細化效果
雙閾值檢測后改進的細化效果
雙閾值檢測后非極大值抑制細化效果