步骤

用高斯滤波器平滑图像；

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

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

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

openCV在C++中的应用

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

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

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

Mat

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

Mat作为一个类，包含

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

从文件中加载图像：

1	Mat img = imread(filename);

如果需要加载灰度图：

1	Mat img = imread(filename, IMREAD_GRAYSCALE);

显示图像：

1 2	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

使用高斯滤镜模糊图像。

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

参量

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

高斯滤波器的C++实现

对图像使用一维高斯卷积模版，在一个方向上进行滤波（例如水平方向）；
转置图像；
对转置后的图像使用同一个高斯卷积模版，在同样的方向上进行滤波；
将图像转置回原来的位置，得到二维高斯滤波的图像。

一维高斯卷积模版可以由二项式展开的系数来模拟，如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)。

a00	a01	a02
a10	a11	a12
a20	a21	a22

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); //计算梯度方向
        }
    }
}

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

三维梯度算法

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

仅仅得到全局梯度并不足以确定边缘，保留局部梯度最大的点，而抑制非极大点。将梯度角的变化范围减小到圆周的四个扇区之一；

四个扇区的标号为0到3，对应3*3领域的四种可能组合方向；
每一个点上领域的中心像素M与沿着梯度线的两个像素比较；
如果M梯度值不比沿梯度线的两个相邻像素梯度值大，则令M=0。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
/**
 局部非极大值抑制
 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

【2】OpenCV Tutorials

【3】OpenCV教程