基于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操作的非极大值抑制有劣势。
双阈值检测后
效果展示
细化前
改进的细化效果
非极大值抑制细化效果
双阈值检测后改进的细化效果
双阈值检测后非极大值抑制细化效果