Canny算子实现边缘特征提取

利用C++与OpenCV库实现Canny算子边缘特征提取，虽然利用OpenCV可以直接调用现有的函数实现边缘提取，但学习嘛，还是要知道一些基础的原理，所以我决定学习学习。。。

文章目录

Canny算子实现边缘特征提取
- 一、基础原理
- 二、检测结果与原码
- 三、结果分析

一、基础原理

Canny算子求边缘点具体步骤可以分为四步：
（如果处理的图像为彩色，那么首先转为灰度图
Gray=0.299R+0.587G+0.114B）
1.用高斯滤波器平滑图像
2.用一阶偏导有限差分计算梯度幅值和放向
3.对梯度幅值进行非极大抑制
4.用双阈值算法检测和连接边缘。

在进行边缘点确定的过程中，并不总是一刀切，也就是超过阈值的点都是边缘点，而是通过设置两个阈值，希望在高阈值和低阈值之间的点也可能是边缘点，而且这些点最好在高阈值附近，也就是说中间阈值的点是高阈值边缘的一种延伸。
对图像进行高斯滤波的实现可以用两个以为高斯核分别两次加权实现，也就是先一维X方向卷积，得到的结果在一维Y方向卷积。也可以直接通过二维高斯核一次卷积完成。
计算梯度值和方向，图像的灰度值的梯度一般使用一阶有限差分来近似，就可以得到图像在X和Y方向上的偏导数的两个矩阵。非极大抑制就是寻找像素点局部最大值，沿着梯度方向，比较它前面和后面的梯度值就行。
双阈值的选取是按照直方图来选择的，首先把梯度幅值的直方图求出来，选取占直方图总数多少所对应的的梯度幅值为高阈值，高阈值的一般为低阈值，也可以选择其他策略。边缘检测时首先判断该点是否超过高阈值，然后判断该点的8邻域点中寻找满足超过低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合。

二、检测结果与原码

在同一实验数据不变的情况下，通过设置多组双阈值（高阈值和低阈值），对实验数据进行边缘特征提取，预期得到在高阈值不变的情况下，通过改变低阈值找到检测效果较好的低阈值。

先放出检测结果：（同时展示出中间结果图）
一共展示出7张图片，分别为：原图、灰度图、平滑后的图像、sobelx、sobley、非极大抑制后的图像、滞后阈值处理后的图像（最终结果）

#include <iostream>
#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <vector>#include "opencv2/imgproc/types_c.h"
#include "opencv2/imgproc/imgproc_c.h"
#include<opencv2/highgui/highgui.hpp>#include <opencv2/highgui/highgui_c.h>
#include<stdio.h>
#include<stdlib.h>
#include<cmath>using namespace std;
using namespace cv;
#define NameWindow "图片处理效果"  //打开的图片路径
string ImagePath = "001.bmp";//卷积运算
void conv2D(InputArray src, InputArray kernel, OutputArray dst, int ddepth, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{//卷积运算的第一步：将卷积核逆时针旋转180度Mat kernelFlip;flip(kernel, kernelFlip, -1);//卷积运算的第二步filter2D(src, dst, ddepth, kernelFlip, anchor, 0.0, borderType);//src输入矩阵，dst输出矩阵，ddepth输出矩阵的数据位深，kernel卷积核（且数据类型为CV_32F/CV_64F），anchor锚点的位置，delta默认值为0，borderType边界扩充的类型。
}
//可分离的same卷积
//可分离的离散二维卷积，先进行垂直方向的卷积，然后进行水平方向上的卷积
void sepConv2D_Y_X(InputArray src, OutputArray src_kerY_ker_X, int ddepth, InputArray kernelY,InputArray kernelX, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{//输入矩阵与垂直方向上的卷积核的卷积Mat src_kerY;conv2D(src, kernelY, src_kerY, ddepth, anchor, borderType);//上面得到到的卷积结果接着和水平方向上的卷积核卷积conv2D(src_kerY, kernelX, src_kerY_ker_X, ddepth, anchor, borderType);//src_kerY_ker_X即为得到的平滑后的结果
}
//先进行水平方向上的卷积，然后进行垂直方向上的卷积
void sepConv2D_X_Y(InputArray src, OutputArray src_kerX_kerY, int ddepth, InputArray kernelX, InputArray kernelY, Point anchor = Point(-1, -1),int borderType = BORDER_DEFAULT)
{//输入矩阵与水平方向上的卷积核的卷积Mat src_kerX;conv2D(src, kernelX, src_kerX, ddepth, anchor, borderType);//上面得到的卷积结果，接着和垂直方向上的卷积核卷积，得到最终的输出结果conv2D(src_kerX, kernelY, src_kerX_kerY, ddepth, anchor, borderType);
}//高斯平滑
Mat gaussBlur(const Mat& image, Size winSize, float sigma, int ddepth = CV_64F, Point anchor = Point(-1, -1), int borderType = BORDER_DEFAULT)
{//卷积核的宽、高均为奇数CV_Assert(winSize.width % 2 == 1 && winSize.height % 2 == 1);//构建垂直方向上的高斯卷积算子Mat gK_y = getGaussianKernel(sigma, winSize.height, CV_64F);//构建水平方向上的高斯卷积算子Mat gK_x = getGaussianKernel(sigma, winSize.width, CV_64F);gK_x = gK_x.t();//转置//分离的高斯卷积Mat blurImage;sepConv2D_Y_X(image, blurImage, ddepth, gK_y, gK_x, Point(-1, -1));return blurImage;
}
//实现非极大抑制
Mat non_maximum_suppressiom_default(Mat dx, Mat dy)
{//使用平方和开方的方式计算边缘强度Mat edgeMag = Mat::zeros(dx.rows, dx.cols, CV_64FC1);magnitude(dx, dy, edgeMag);//得到的edgeMag即为计算好的边缘强度矩阵，此函数为OpenCV//下面用自己写出来的计算边缘强度的效果没有magnitude()函数的效果好，不过差距不大，但还是选择使用magnitude()函数计算/*for (int i = 0;i < dx.rows; i++){for (int j = 0;j < dx.cols;j++){edgeMag.at<double>(i, j) = sqrt(pow(dx.at<double>(i, j), 2) + pow(dy.at<double>(i, i), 2));}}*///高、宽Mat edgeMag_nonMaxSup = Mat::zeros(dx.size(), dx.type());//构造一个与dx行列数相同、数据类型相同的全是0的矩阵int rows = dx.rows;//这里选用dx、dy都可以int cols = dx.cols;//边缘强度的非极大抑制//如果magnitude(r,c)在沿着梯度方向angle(r,c)上的邻域内是最大的则为极大值；否则，设置为0。//Mat angle = Mat::zeros(dx.rows, dx.cols, CV_8UC1);for (int r = 1; r < rows - 1; r++){for (int c = 1; c < cols - 1; c++){double x = dx.at<double>(r, c);double y = dy.at<double>(r, c);//梯度方向double angle = atan2(y, x / CV_PI * 180);//当前位置的边缘强度float mag = edgeMag.at<double>(r, c);//左、右方向比较if (abs(angle) < 22.5 || abs(angle) > 157.5){double left = edgeMag.at<double>(r, c - 1);double right = edgeMag.at<double>(r, c + 1);if (mag > left && mag > right)edgeMag_nonMaxSup.at<double>(r, c) = mag;}//左上、右下方向比较if ((angle >= 22.5 && angle < 67.5) || (angle < -112.5 && angle >= -157.5)){double leftTop = edgeMag.at<double>(r - 1, c - 1);double rightBottom = edgeMag.at<double>(r + 1, c + 1);if (mag > leftTop && mag > rightBottom){edgeMag_nonMaxSup.at<double>(r, c) = mag;}}//上、下方向比较if ((angle >= 67.5 && angle <= 112.5) || (angle >= -112.5 && angle <= -67.5)){double top = edgeMag.at<double>(r - 1, c);double bottom = edgeMag.at<double>(r + 1, c);if (mag > top && mag > bottom){edgeMag_nonMaxSup.at<double>(r, c) = mag;}}//右上、左下方向比较if ((angle > 112.5 && angle <= 157.5) || (angle >= -67.5 && angle <= -22.5)){double rightTop = edgeMag.at<double>(r - 1, c + 1);double leftBottom = edgeMag.at<double>(r + 1, c - 1);if (mag > rightTop && mag > leftBottom){edgeMag_nonMaxSup.at<double>(r, c) = mag;}}}}return edgeMag_nonMaxSup;
}
//滞后阈值处理
//确定一个点的坐标是否在图像范围内
bool checkInRange(int r, int c, int rows, int cols)
{if (r >= 0 && r < rows && c >= 0 && c <= cols)return true;elsereturn false;
}
//从确定边缘点出发，延长边缘
//void trace(Mat edgeMag_nonMaxSup, Mat& edge, float lowerThresh, int r, int c, int rows, int cols)
//{//  if (edge.at<uchar>(r, c) == 0)
//  {//      edge.at<uchar>(r, c) = 255;
//      for (int i = -1;i <= 1;i++)
//      {//          for (int j = -1;j <= 1;j++)
//          {//              double mag = edgeMag_nonMaxSup.at<double>(r + i, c + i);
//              if (checkInRange(r + i, c + j, rows, cols) && mag >= lowerThresh)
//                  trace(edgeMag_nonMaxSup, edge, lowerThresh, r + i, c + j,rows, cols);
//          }
//      }
//  }
//}
//双阈值的滞后阈值处理
//Mat hysteresisThreshold(Mat edgeMag_nonMaxSup, float lowerThresh, float upperThresh)
//{//  //高、宽
//  int rows = edgeMag_nonMaxSup.rows;
//  int cols = edgeMag_nonMaxSup.cols;
//  //最后的边缘输出图
//  Mat edge = Mat::zeros(Size(cols, rows), CV_64FC1);
//  //滞后阈值处理
//  for (int r = 1;r < rows - 1;r++)
//  {//      for (int c = 1;r < cols - 1;c++)
//      {//          double mag=edgeMag_nonMaxSup.at<double>(r, c);
//          //大于阈值的点，可作为确定边缘点被接受
//          //并以该点为起始点延长边缘
//          if (mag >= upperThresh)
//              //edge.convertTo(edge, CV_8UC3);
//              trace(edgeMag_nonMaxSup, edge, lowerThresh, r, c, rows, cols);
//          //小于低阈值的点直接被剔除
//          if (mag < lowerThresh)
//              edge.at<double>(r, c) = 0;
//      }
//  }
//  return edge;
//}Mat hysteresisThresholdvalue(Mat edgeMag_nonMaxSup, float lowerThresh, float upperThresh)
{//高、宽int rows = edgeMag_nonMaxSup.rows;int cols = edgeMag_nonMaxSup.cols;Mat edge = Mat::zeros(edgeMag_nonMaxSup.size(), CV_64FC1);for (int r = 1;r < rows - 1;r++){for (int c = 1;c < cols - 1;c++){double mag = edgeMag_nonMaxSup.at<double>(r, c);//在（r，c）点的灰度值if (mag < lowerThresh) //当该点的灰度值小于最小阈值时直接剔除edge.at<double>(r, c) = 0;if (mag > upperThresh)//当该点的灰度值大于最大阈值时确定为边缘点edge.at<double>(r, c) = 255;if (mag <= upperThresh && mag >= lowerThresh)//当该点的灰度值在最大最小阈值之间时，则判断其8领域内for (int i = -1;i <= 1;i++)              //是否有一个点的灰度值大于最小阈值，若有则该点也判断为边缘点{for (int j = -1;j <= 1;j++){double mag1 = edgeMag_nonMaxSup.at<double>(r + i, c + i);if (checkInRange(r + i, c + j, rows, cols) && (mag1 >= mag))edge.at<double>(r, c) = 255;}}}}return edge;
}//自定义一个窗口显示多图函数
void imshowMany(const std::string& _winName, const vector<Mat>& ployImages)
{int nImg = (int)ployImages.size();//获取在同一画布中显示多图的数目  Mat dispImg;int size;int x, y;//w：多图按矩阵排列的行数  ，h: 多图按矩阵排列的的列数    int w, h;float scale;//缩放比例  int max;if (nImg <= 0){printf("Number of arguments too small....\n");return;}else if (nImg > 12){printf("Number of arguments too large....\n");return;}else if (nImg == 1){w = h = 1;size = 400;}else if (nImg == 2){w = 2; h = 1;//2x1  size = 400;}else if (nImg == 3 || nImg == 4){w = 2; h = 2;//2x2  size = 450;}else if (nImg == 5 || nImg == 6){w = 3; h = 2;//3x2  size = 300;}else if (nImg == 7 || nImg == 8){w = 4; h = 2;//4x2  size = 300;}else{w = 4; h = 3;//4x3  size = 200;}dispImg.create(Size(100 + size * w, 30 + size * h), CV_8UC3);//根据图片矩阵w*h，创建画布，可线的图片数量为w*h  for (int i = 0, m = 20, n = 20; i < nImg; i++, m += (20 + size)){x = ployImages[i].cols; //第(i+1)张子图像的宽度(列数)  y = ployImages[i].rows;//第(i+1)张子图像的高度（行数）  max = (x > y) ? x : y;//比较每张图片的行数和列数，取大值  scale = (float)((float)max / size);//计算缩放比例  if (i % w == 0 && m != 20){m = 20;n += 20 + size;}Mat imgROI = dispImg(Rect(m, n, (int)(x / scale), (int)(y / scale))); //在画布dispImage中划分ROI区域  resize(ployImages[i], imgROI, Size((int)(x / scale), (int)(y / scale))); //将要显示的图像设置为ROI区域大小  }namedWindow(_winName);imshow(_winName, dispImg);
}int main()
{ Mat image = imread("888.jpg");image.convertTo(image, CV_64FC3);/*图片默认是unsigned char型，范围是0-255*/Mat image_gray = Mat::zeros(image.rows, image.cols, CV_64FC1);Mat image_smooth = Mat::zeros(image.rows, image.cols, CV_64FC1);Mat dx = Mat::zeros(image.rows, image.cols, CV_64FC1);Mat dy = Mat::zeros(image.rows, image.cols, CV_64FC1);Mat image_restrain = Mat::zeros(image.rows, image.cols, CV_64FC1);Mat dst = Mat::zeros(image.rows, image.cols, CV_64FC1);//水平方向卷积核sobelxMat sobelx = (Mat_<double>(3, 3) << 1, 0, -1,2, 0, -2,1, 0, -1);//垂直方向卷积核sobelyMat  sobely = (Mat_<double>(3, 3) << 1, 2, 1,0, 0, 0,-1, -2, -1);/*++++++++++++++++++++++++++++++++++++++++++++图像灰度化++++++++++++++++++++++++++++++++++*/for (int i = 0; i < image.rows; i++) {for (int j = 0; j < image.cols; j++) {image_gray.at<double>(i, j) = (0.114 * image.at<Vec3d>(i, j)[0] + 0.299 * image.at<Vec3d>(i, j)[1] + 0.587 * image.at<Vec3d>(i, j)[2]);}}image.convertTo(image, CV_8UC3);/*+++++++++++++++++++++++++++++++++++++高斯滤波平滑处理++++++++++++++++++++++++++++++++++++++++++++*/image_smooth=gaussBlur(image_gray, Size(5,5), 3, -1,Point(-1,-1),BORDER_CONSTANT);//GaussianBlur(image_gray, image_smooth, Size(5, 5), 3, 0, BORDER_DEFAULT);//OpenCV实现高斯平滑的函数/*++++++++++++++++++++++++++++++++++++++++非极大抑制++++++++++++++++++++++++++++++++++++++++++++++++*/conv2D(image_smooth, sobelx, dx, CV_64FC1, Point(-1, -1), BORDER_CONSTANT);conv2D(image_smooth, sobely, dy, CV_64FC1, Point(-1, -1), BORDER_CONSTANT);image_restrain=non_maximum_suppressiom_default(dx, dy);/*++++++++++++++++++++++++++++++++++++++滞后阈值处理++++++++++++++++++++++++++++++++++++++++++++++++++*/int min, max;cout << "输入阈值范围:" << endl;cin >> min >> max;dst = hysteresisThresholdvalue(image_restrain, min, max);image_gray.convertTo(image_gray, CV_8UC1);Mat Image2,Image3,Image4,Image5,Image6,Image7;cvtColor(image_gray, Image2, CV_GRAY2BGR);//转为三通道图,用于多图显示，（因为只能显示图），三通道位置（1，2）image_smooth.convertTo(image_smooth, CV_8UC1);cvtColor(image_smooth, Image3, CV_GRAY2BGR);dx.convertTo(dx, CV_8UC1);cvtColor(dx, Image4, CV_GRAY2BGR);dy.convertTo(dy, CV_8UC1);cvtColor(dy, Image5, CV_GRAY2BGR);image_restrain.convertTo(image_restrain, CV_8UC1);cvtColor(image_restrain, Image6, CV_GRAY2BGR);dst.convertTo(dst, CV_8UC1);cvtColor(dst, Image7, CV_GRAY2BGR);vector<Mat>Images(7);//模板类vector，用于放置7个类型为Mat的元素，即七张图片  Images[0] = image;Images[1] = Image2;Images[2] = Image3;Images[3] = Image4;Images[4] = Image5;Images[5] = Image6;Images[6] = Image7;namedWindow(NameWindow);imshowMany(NameWindow, Images);//调用一个窗口显示多图函数 /*imshow("原图", image);imshow("灰度图", image_gray);imshow("平滑后图像", image_smooth);imshow("sobelx", dx);imshow("sobley", dy);imshow("非极大抑制", image_restrain);imshow("滞后阈值处理", dst);*/waitKey(0);return 0;
}

三、结果分析

1.设置高阈值为150不变，改变低阈值

2.设置低阈值为30不变，改变高阈值

3.分析

当实验中选取高阈值150固定不变，低阈值分别选取了30、50、75、100对原图像进行了Canny算法的特征检测及提取。通过对不同的检测结果分析，可以看出，对于同一数据，在高阈值不变的情况下，低阈值越低检测结果中噪点越多，检测的非边缘点越多，随着低阈值的变大，检测结果中噪点明显减少，但当低阈值增大到一定程度，对部分边缘的检测就会遗失。经分析当低阈值与高阈值之间的比值为1:3是检测效果最佳。

但当实验中选取低阈值为30固定不变，高阈值分别选取了45、60、75、90对原图像进行了Canny算法的特征检测及提取。通过对不同的检测结果分析，可以看出，变化并不大，所有对于同一数据，高、低阈值不仅要选择最适合的比值，高、低阈值大小的选择也很重要。

双阈值算法检测中，高、阈值要选择最恰当的比例关系，也要根据图像特点，选择合适的阈值范围。