边缘检测算法的C++实现

基于sobel算子的边缘检测算法的C++实现及改进

实验内容

  1. 对已有的sobel算子边缘检测算法的学习
  2. 对已有算法进行C++实现
  3. 对已有算法进行改进

已有的边缘检测算法流程

  1. 对图像进行灰度化得到灰度图像
  2. 对灰度图像进行高斯模糊去除部分噪声
  3. 利用sobel算子作为卷积核对图像进行卷积，计算梯度
  4. 用非极大值抑制法进行边缘细化
  5. 利用双阈值进行噪声去除

算法实现

图像灰度化
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}\)

其中G即为某个像素点的梯度值，而梯度方向则为 $arctan({G_y \over G_x})$ 
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);
}
固定单阈值

然而固定单阈值法有许多问题

  1. 对于不同类型的图像，不同要求的边缘检测，阈值不同，得到效果好的阈值很麻烦
  2. 极容易造成边缘的中断
    因此采取双阈值法可以一定程度上解决问题

我们定义高阈值和低阈值，对于灰度值高于高阈值的像素点我们将其确定为边缘，而灰度值低于低阈值的像素点我们认为它一定不是边缘，对于灰度值介于高阈值和低阈值之间的像素点，它有可能是边缘，当它相邻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操作的非极大值抑制有劣势。

双阈值检测后

效果展示

细化前
改进的细化效果
非极大值抑制细化效果
双阈值检测后改进的细化效果
双阈值检测后非极大值抑制细化效果
posted @ 2019-12-17 15:16  thhyj  阅读(2811)  评论(0)  编辑  收藏  举报