Canny边缘检测算法(C++实现)

步骤

  1. 用高斯滤波器平滑图像;
  2. 用一阶偏导有限差分计算梯度幅值和方向;
  3. 对梯度幅值应用非极大值抑制;
  4. 用双阈值算法检测和连接边缘。

 

一、openCV在C++中的应用

首先,在mac的Xcode上安装配置openCV库,参考一下链接(科学上网访问)https://medium.com/@jaskaranvirdi/setting-up-opencv-and-c-development-environment-in-xcode-b6027728003

using namespace cv; // 使用命名空间cv

如此可以减少输入,例如 cv :: Mat 就可省略为 Mat

Mat

Mat的优点是不再需要手动分配其内存,并在不需要它时立即发布它。在执行此操作仍然是可能的情况下,大多数OpenCV功能将自动分配其输出数据。

Mat作为一个类,包含

  • 矩阵头(包含矩阵的大小,用于存储的方法,存储在哪个地址的信息等等)
  • 指向包含像素值(取决于所选存储方法的任何维度)

从文件中加载图像:

Mat img = imread(filename);

如果需要加载灰度图:

Mat img = imread(filename, IMREAD_GRAYSCALE);

显示图像:

namedWindow("图片"); //打开名为“图片”的窗口
imshow("图片", img); //显示图像

 openCV加载图像显示

using namespace cv; // 使用命名空间cv

int main()
{
    Mat img = imread("building.jpg", IMREAD_GRAYSCALE); //从文件中加载灰度图像
    //显示图像
    namedWindow("原图");
    imshow("原图", img);

    waitKey(); //等待键值输入
    return 0;
}

如何访问图像每一个像素点

利用指针访问,调用 Mat::ptr(i) 来获取第i行的首地址,通过循环进行访问。

// 按行遍历所有点(单通道)
for (int j = 0; j < nr; j++) {
    for (int i = 0; i < nc; i++) {
        //每个点为 img.ptr<uchar>(j)[i]
    }
}

 

二、用高斯滤波器平滑图像

高斯滤波器(openCV)

openCV自带的高斯滤波器:cv :: GaussianBlur

void cv::GaussianBlur (InputArray src
  OutputArray dst
  Size ksize
  double sigmaX
  double sigmaY = 0
  int borderType = BORDER_DEFAULT 
 )

#include <opencv2/imgproc.hpp>

使用高斯滤镜模糊图像。

该函数将源图像与指定的高斯内核进行卷积。支持就地过滤。

参量

src输入图像;图像可以具有任意数量的通道,这些通道可以独立处理,但深度应为CV_8U,CV_16U,CV_16S,CV_32F或CV_64F。
dst输出与src大小和类型相同的图像。
size高斯核大小。ksize.width和ksize.height可以不同,但​​它们都必须为正数和奇数。或者,它们可以为零,然后根据sigma计算得出。
sigmaXX方向上的高斯核标准偏差。
sigmaYY方向的高斯核标准差;如果sigmaY为零,则将其设置为等于sigmaX;如果两个sigmas为零,则分别从ksize.width和ksize.height计算得出(有关详细信息,请参见getGaussianKernel);为了完全控制结果,而不管将来可能对所有这些语义进行的修改,建议指定所有ksize,sigmaX和sigmaY。
borderType像素外推方法,请参见BorderTypes

高斯滤波器的C++实现

  1. 对图像使用一维高斯卷积模版,在一个方向上进行滤波(例如水平方向);
  2. 转置图像;
  3. 对转置后的图像使用同一个高斯卷积模版,在同样的方向上进行滤波;
  4. 将图像转置回原来的位置,得到二维高斯滤波的图像。

一维高斯卷积模版可以由二项式展开的系数来模拟,如3*3模版:     1/4 * [1  2  1]

/**
 高斯滤波器,利用3*3的高斯模版进行高斯卷积
 img 输入原图像
 dst  高斯滤波后的输出图像
*/
void gaussianFilter(Mat &img, Mat &dst) {
    // 对水平方向进行滤波
    Mat dst1 = img.clone();
    gaussianConvolution(img, dst1);
    //图像矩阵转置
    Mat dst2;
    transpose(dst1, dst2);
    // 对垂直方向进行滤波
    Mat dst3 = dst2.clone();
    gaussianConvolution(dst2, dst3);
    // 再次转置
    transpose(dst3, dst);
}
/**
 一维高斯卷积,对每行进行高斯卷积
 img 输入原图像
 dst  一维高斯卷积后的输出图像
 */
void gaussianConvolution(Mat &img, Mat &dst) {
    int nr = img.rows;
    int nc = img.cols;
    int templates[3] = {1, 2, 1};
    
    // 按行遍历除每行边缘点的所有点
    for (int j = 0; j < nr; j++) {
        uchar* data= img.ptr<uchar>(j); //提取该行地址
        for (int i = 1; i < nc-1; i++) {
            int sum = 0;
            for (int n = 0; n < 3; n++) {
                sum += data[i-1+n] * templates[n]; //相称累加
            }
            sum /= 4;
            dst.ptr<uchar>(j)[i] = sum;
        }
    }
}

高斯滤波前后的图像:

高斯滤波

 

 三、用一阶偏导有限差分计算梯度幅值和方向

用一阶偏导有限差分计算偏导数的两个阵列P与Q

P\left [ y,x \right ] \approx \left ( S\left [ y,x+1 \right ]-S\left [ y,x \right ] +S\left [ y+1,x+1 \right ]-S[y+1,x]\right )/2

Q\left [ y,x \right ] \approx \left ( S\left [ y+1,x \right ]-S\left [ y,x \right ] +S\left [ y+1,x+1 \right ]-S[y,x+1]\right )/2

再由P和Q算出梯度幅值和方向角

M\left [ y,x \right ]= \sqrt{P\left [ y,x \right ]^{2}+Q\left [ y,x \right ]^{2}}

\theta \left [ y,x \right ]=\arctan \left ( Q\left [ y,x \right ] /P\left [ y,x \right ]\right )

/**
 用一阶偏导有限差分计算梯度幅值和方向
 img 输入原图像
 gradXY 输出的梯度幅值
 theta 输出的梯度方向
 */
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
    gradXY = Mat::zeros(img.size(), CV_8U);
    theta = Mat::zeros(img.size(), CV_8U);
    
    for (int i = 0; i < img.rows-1; i++) {
        for (int j = 0; j < img.cols-1; j++) {
            double p = (img.ptr<uchar>(j)[i+1] - img.ptr<uchar>(j)[i] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j+1)[i])/2;
            double q = (img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j)[i] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j)[i+1])/2;
            gradXY.ptr<uchar>(j)[i] = sqrt(p*p + q*q); //计算梯度
            theta.ptr<uchar>(j)[i] = atan(q/p);
        }
    }
}

此时输入输出图像为:

二维梯度算法

 可以看出,二维计算梯度只区分出了部分边界,边界损失过大,于是采用三维算法计算梯度((y,x)为a11)。

a00a01a02
a10a11a12
a20a21a22
double gradX = double(a02 + 2 * a12 + a22 - a00 - 2 * a10 - a20);
double gradY = double(a00 + 2 * a01 + a02 - a20 - 2 * a21 - a22);
/**
 用一阶偏导有限差分计算梯度幅值和方向
 img 输入原图像
 gradXY 输出的梯度幅值
 theta 输出的梯度方向
 */
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
    gradXY = Mat::zeros(img.size(), CV_8U);
    theta = Mat::zeros(img.size(), CV_8U);
    
    for (int j = 1; j < img.rows-1; j++) {
        for (int i = 1; i < img.cols-1; i++) {
            double gradY = double(img.ptr<uchar>(j-1)[i-1] + 2 * img.ptr<uchar>(j-1)[i] + img.ptr<uchar>(j-1)[i+1] - img.ptr<uchar>(j+1)[i-1] - 2 * img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j+1)[i+1]);
            double gradX = double(img.ptr<uchar>(j-1)[i+1] + 2 * img.ptr<uchar>(j)[i+1] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j-1)[i-1] - 2 * img.ptr<uchar>(j)[i-1] - img.ptr<uchar>(j+1)[i-1]);

            gradXY.ptr<uchar>(j)[i] = sqrt(gradX*gradX + gradY*gradY); //计算梯度
            theta.ptr<uchar>(j)[i] = atan(gradY/gradX); //计算梯度方向
        }
    }
}

三维梯度算法的输入输出图像:

三维梯度算法

 

四、对梯度幅值应用非极大值抑制

仅仅得到全局梯度并不足以确定边缘,保留局部梯度最大的点,而抑制非极大点。

  1. 将梯度角的变化范围减小到圆周的四个扇区之一;
  2. 四个扇区的标号为0到3,对应3*3领域的四种可能组合方向;
  3. 每一个点上领域的中心像素M与沿着梯度线的两个像素比较;
  4. 如果M梯度值不比沿梯度线的两个相邻像素梯度值大,则令M=0。    

      由 atan() 得到的角度在 \left [ -\pi /2, \pi /2\right ] 范围内,将此范围均分为四个等份。

 

                                                           

/**
 局部非极大值抑制
 gradXY 输入的梯度幅值
 theta 输入的梯度方向
 dst 输出的经局部非极大值抑制后的图像
 */
void nonLocalMaxValue (Mat &gradXY, Mat &theta, Mat &dst) {
    dst = gradXY.clone();
    for (int j = 1; j < gradXY.rows-1; j++) {
        for (int i = 1; i < gradXY.cols-1; i++) {
            double t = double(theta.ptr<uchar>(j)[i]);
            double g = double(dst.ptr<uchar>(j)[i]);
            if (g == 0.0) {
                continue;
            }
            double g0, g1;
            if ((t >= -(3*M_PI/8)) && (t < -(M_PI/8))) {
                g0 = double(dst.ptr<uchar>(j-1)[i-1]);
                g1 = double(dst.ptr<uchar>(j+1)[i+1]);
            }
            else if ((t >= -(M_PI/8)) && (t < M_PI/8)) {
                g0 = double(dst.ptr<uchar>(j)[i-1]);
                g1 = double(dst.ptr<uchar>(j)[i+1]);
            }
            else if ((t >= M_PI/8) && (t < 3*M_PI/8)) {
                g0 = double(dst.ptr<uchar>(j-1)[i+1]);
                g1 = double(dst.ptr<uchar>(j+1)[i-1]);
            }
            else {
                g0 = double(dst.ptr<uchar>(j-1)[i]);
                g1 = double(dst.ptr<uchar>(j+1)[i]);
            }
            
            if (g <= g0 || g <= g1) {
                dst.ptr<uchar>(j)[i] = 0.0;
            }
        }
    }
}

输入的经梯度计算后的图像和输出的局部非极大值抑制后的图像:

局部非极大值抑制

五、用双阈值算法检测和连接边缘

1、Canny算法采用双阈值,高阈值一般是低阈值的两倍,遍历所有像素点:

X < 低阈值 ,像素点置0,被抑制掉;

低阈值 < X <高阈值,像素点为弱边缘点,像素点值先不变;

X > 高阈值,像素点为强边缘点,置255。

2、弱边缘点补充连接强边缘点:

如果弱边缘点的8邻点域存在强边缘点,则将此点置255,用以连接强边缘点;如果不存在强边缘点,则这是一个孤立的弱边缘点,此点置0。

/**
 用双阈值算法检测和连接边缘
 low 输入的低阈值
 high 输入的高阈值
 img 输入的原图像
 dst 输出的用双阈值算法检测和连接边缘后的图像
 */
void doubleThreshold (double low, double high, Mat &img, Mat &dst) {
    dst = img.clone();
    
    // 区分出弱边缘点和强边缘点
    for (int j = 0; j < img.rows-1; j++) {
        for (int i = 0; i < img.cols-1; i++) {
            double x = double(dst.ptr<uchar>(j)[i]);
            // 像素点为强边缘点,置255
            if (x > high) {
                dst.ptr<uchar>(j)[i] = 255;
            }
            // 像素点置0,被抑制掉
            else if (x < low) {
                dst.ptr<uchar>(j)[i] = 0;
            }
        }
    }
    
    // 弱边缘点补充连接强边缘点
    doubleThresholdLink(dst);
}
/**
 弱边缘点补充连接强边缘点
 img 弱边缘点补充连接强边缘点的输入和输出图像
 */
void doubleThresholdLink (Mat &img) {
    // 循环找到强边缘点,把其领域内的弱边缘点变为强边缘点
    for (int j = 1; j < img.rows-2; j++) {
        for (int i = 1; i < img.cols-2; i++) {
            // 如果该点是强边缘点
            if (img.ptr<uchar>(j)[i] == 255) {
                // 遍历该强边缘点领域
                for (int m = -1; m < 1; m++) {
                    for (int n = -1; n < 1; n++) {
                        // 该点为弱边缘点(不是强边缘点,也不是被抑制的0点)
                        if (img.ptr<uchar>(j+m)[i+n] != 0 && img.ptr<uchar>(j+m)[i+n] != 255) {
                            img.ptr<uchar>(j+m)[i+n] = 255; //该弱边缘点补充为强边缘点
                        }
                    }
                }
            }
        }
    }
    
    for (int j = 0; j < img.rows-1; j++) {
        for (int i = 0; i < img.cols-1; i++) {
            // 如果该点依旧是弱边缘点,及此点是孤立边缘点
            if (img.ptr<uchar>(j)[i] != 255 && img.ptr<uchar>(j)[i] != 255) {
                img.ptr<uchar>(j)[i] = 0; //该孤立弱边缘点抑制
            }
        }
    }
}

双阈值算法前后的输入输出图像 :

双阈值算法

 

Canny边缘检测代码

#include <opencv2/opencv.hpp>
#include <math.h>

#define _USE_MATH_DEFINES

using namespace cv;// 使用命名空间cv


/**
 将两个图像拼接,以便在同一个窗口显示
 dst 输出的拼接后的图像
 src1 拼接的第一幅图
 src2 拼接的第二幅图
 */
void mergeImg(Mat & dst,Mat &src1,Mat &src2) {
    int rows = src1.rows;
    int cols = src1.cols+5+src2.cols;
    CV_Assert(src1.type () == src2.type ());
    dst.create (rows,cols,src1.type ());
    src1.copyTo (dst(Rect(0,0,src1.cols,src1.rows)));
    src2.copyTo (dst(Rect(src1.cols+5,0,src2.cols,src2.rows)));
}


/**
 一维高斯卷积,对每行进行高斯卷积
 img 输入原图像
 dst  一维高斯卷积后的输出图像
 */
void gaussianConvolution(Mat &img, Mat &dst) {
    int nr = img.rows;
    int nc = img.cols;
    int templates[3] = {1, 2, 1};
    
    // 按行遍历除每行边缘点的所有点
    for (int j = 0; j < nr; j++) {
        uchar* data= img.ptr<uchar>(j); //提取该行地址
        for (int i = 1; i < nc-1; i++) {
            int sum = 0;
            for (int n = 0; n < 3; n++) {
                sum += data[i-1+n] * templates[n]; //相称累加
            }
            sum /= 4;
            dst.ptr<uchar>(j)[i] = sum;
        }
    }
}


/**
 高斯滤波器,利用3*3的高斯模版进行高斯卷积
 img 输入原图像
 dst  高斯滤波后的输出图像
*/
void gaussianFilter(Mat &img, Mat &dst) {
    // 对水平方向进行滤波
    Mat dst1 = img.clone();
    gaussianConvolution(img, dst1);
    //图像矩阵转置
    Mat dst2;
    transpose(dst1, dst2);
    // 对垂直方向进行滤波
    Mat dst3 = dst2.clone();
    gaussianConvolution(dst2, dst3);
    // 再次转置
    transpose(dst3, dst);
}


/**
 用一阶偏导有限差分计算梯度幅值和方向
 img 输入原图像
 gradXY 输出的梯度幅值
 theta 输出的梯度方向
 */
void getGrandient (Mat &img, Mat &gradXY, Mat &theta) {
    gradXY = Mat::zeros(img.size(), CV_8U);
    theta = Mat::zeros(img.size(), CV_8U);
    
    for (int j = 1; j < img.rows-1; j++) {
        for (int i = 1; i < img.cols-1; i++) {
            double gradY = double(img.ptr<uchar>(j-1)[i-1] + 2 * img.ptr<uchar>(j-1)[i] + img.ptr<uchar>(j-1)[i+1] - img.ptr<uchar>(j+1)[i-1] - 2 * img.ptr<uchar>(j+1)[i] - img.ptr<uchar>(j+1)[i+1]);
            double gradX = double(img.ptr<uchar>(j-1)[i+1] + 2 * img.ptr<uchar>(j)[i+1] + img.ptr<uchar>(j+1)[i+1] - img.ptr<uchar>(j-1)[i-1] - 2 * img.ptr<uchar>(j)[i-1] - img.ptr<uchar>(j+1)[i-1]);

            gradXY.ptr<uchar>(j)[i] = sqrt(gradX*gradX + gradY*gradY); //计算梯度
            theta.ptr<uchar>(j)[i] = atan(gradY/gradX); //计算梯度方向
        }
    }
}


/**
 局部非极大值抑制
 gradXY 输入的梯度幅值
 theta 输入的梯度方向
 dst 输出的经局部非极大值抑制后的图像
 */
void nonLocalMaxValue (Mat &gradXY, Mat &theta, Mat &dst) {
    dst = gradXY.clone();
    for (int j = 1; j < gradXY.rows-1; j++) {
        for (int i = 1; i < gradXY.cols-1; i++) {
            double t = double(theta.ptr<uchar>(j)[i]);
            double g = double(dst.ptr<uchar>(j)[i]);
            if (g == 0.0) {
                continue;
            }
            double g0, g1;
            if ((t >= -(3*M_PI/8)) && (t < -(M_PI/8))) {
                g0 = double(dst.ptr<uchar>(j-1)[i-1]);
                g1 = double(dst.ptr<uchar>(j+1)[i+1]);
            }
            else if ((t >= -(M_PI/8)) && (t < M_PI/8)) {
                g0 = double(dst.ptr<uchar>(j)[i-1]);
                g1 = double(dst.ptr<uchar>(j)[i+1]);
            }
            else if ((t >= M_PI/8) && (t < 3*M_PI/8)) {
                g0 = double(dst.ptr<uchar>(j-1)[i+1]);
                g1 = double(dst.ptr<uchar>(j+1)[i-1]);
            }
            else {
                g0 = double(dst.ptr<uchar>(j-1)[i]);
                g1 = double(dst.ptr<uchar>(j+1)[i]);
            }
            
            if (g <= g0 || g <= g1) {
                dst.ptr<uchar>(j)[i] = 0.0;
            }
        }
    }
}


/**
 弱边缘点补充连接强边缘点
 img 弱边缘点补充连接强边缘点的输入和输出图像
 */
void doubleThresholdLink (Mat &img) {
    // 循环找到强边缘点,把其领域内的弱边缘点变为强边缘点
    for (int j = 1; j < img.rows-2; j++) {
        for (int i = 1; i < img.cols-2; i++) {
            // 如果该点是强边缘点
            if (img.ptr<uchar>(j)[i] == 255) {
                // 遍历该强边缘点领域
                for (int m = -1; m < 1; m++) {
                    for (int n = -1; n < 1; n++) {
                        // 该点为弱边缘点(不是强边缘点,也不是被抑制的0点)
                        if (img.ptr<uchar>(j+m)[i+n] != 0 && img.ptr<uchar>(j+m)[i+n] != 255) {
                            img.ptr<uchar>(j+m)[i+n] = 255; //该弱边缘点补充为强边缘点
                        }
                    }
                }
            }
        }
    }
    
    for (int j = 0; j < img.rows-1; j++) {
        for (int i = 0; i < img.cols-1; i++) {
            // 如果该点依旧是弱边缘点,及此点是孤立边缘点
            if (img.ptr<uchar>(j)[i] != 255 && img.ptr<uchar>(j)[i] != 255) {
                img.ptr<uchar>(j)[i] = 0; //该孤立弱边缘点抑制
            }
        }
    }
}


/**
 用双阈值算法检测和连接边缘
 low 输入的低阈值
 high 输入的高阈值
 img 输入的原图像
 dst 输出的用双阈值算法检测和连接边缘后的图像
 */
void doubleThreshold (double low, double high, Mat &img, Mat &dst) {
    dst = img.clone();
    
    // 区分出弱边缘点和强边缘点
    for (int j = 0; j < img.rows-1; j++) {
        for (int i = 0; i < img.cols-1; i++) {
            double x = double(dst.ptr<uchar>(j)[i]);
            // 像素点为强边缘点,置255
            if (x > high) {
                dst.ptr<uchar>(j)[i] = 255;
            }
            // 像素点置0,被抑制掉
            else if (x < low) {
                dst.ptr<uchar>(j)[i] = 0;
            }
        }
    }
    
    // 弱边缘点补充连接强边缘点
    doubleThresholdLink(dst);
}


int main () {
    Mat img = imread("woman.jpg", IMREAD_GRAYSCALE); //从文件中加载灰度图像

    // 读取图片失败,则停止
    if (img.empty()) {
        printf("读取图像文件失败");
        system("pause");
        return 0;
    }
    
    // 高斯滤波
    Mat gauss_img;
    gaussianFilter(img, gauss_img); //高斯滤波器
    
    // 用一阶偏导有限差分计算梯度幅值和方向
    Mat gradXY, theta;
    getGrandient(gauss_img, gradXY, theta);
    
    // 局部非极大值抑制
    Mat local_img;
    nonLocalMaxValue(gradXY, theta, local_img);
    
    // 用双阈值算法检测和连接边缘
    Mat dst;
    doubleThreshold(40, 80, local_img, dst);

    // 图像显示
    Mat outImg;
    mergeImg (outImg,img,dst); //图像拼接
    namedWindow("img");
    imshow("img",outImg);// 图像显示
    imwrite("canny算法.jpg", outImg);

    waitKey(); //等待键值输入
    return 0;
}

Canny边缘检测的前后图像

Canny边缘检测算法

 

参考资源:

【1】Setting up OpenCV and C++ development environment in Xcode for Computer Vision projects,地址:https://medium.com/@jaskaranvirdi/setting-up-opencv-and-c-development-environment-in-xcode-b6027728003

【2】OpenCV Tutorials,地址:https://docs.opencv.org/master/d9/df8/tutorial_root.html

【3】OpenCV教程,地址:https://www.w3cschool.cn/opencv/opencv-2gnx28u3.html

 

 

  • 18
    点赞
  • 133
    收藏
    觉得还不错? 一键收藏
  • 3
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值