自适应阈值canny边缘检测（功能实现）

学习记录…

1 概述

canny边缘检测是一种特别常用且性能优秀的边缘检测算法，相比于普通的边缘检测算法，canny获得的边缘较细且具有连续的边缘轮廓，为之后的一系列图像处理带来极大的便利。

canny边缘检测也是基于梯度图像的，通常在其局部最大值附近会包含一些宽脊，为了细化这些宽脊采用的方向就是非极大值抑制——梯度的本意是一个向量（矢量），函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模——即梯度图像像素值），梯度的方向是与边缘的方向垂直的，那么在一个3x3范围内，可以将梯度的方向进行分区：

对于每个像素点，如果 A ( i , j ) A(i,j) A(i,j)的梯度幅值比其梯度方向上相邻2个像素点 A 1 A1 A1和 A 2 A2 A2的梯度幅值大，该点标记为候选边缘点。

梯度方向（角度）在不同的分区可以分别映射为水平方向（垂直边缘）、+45方向、垂直方向（水平边缘）、-45方向。

那么在确定某一点梯度方向所属分区所映射到的方向之后，就将该点梯度幅值与方向上的梯度幅值进行比较，若该点梯度幅值均大于方向上点的梯度幅值则保留，否则令为0。

改进

在canny边缘检测中，还有一个重要的步骤：双阈值的滞后阈值处理，一个高阈值TH和一个低阈值TL，比例在2:1到3:1内，（至于为什么会这样真不明白）这就带来了canny边缘检测的一个很大的缺点，那就是需要输入阈值参数，基于此，很多完全自适应阈值的canny算法诞生，在这里仅提供一种较简单和实用的思路——将经过非极大值抑制后的梯度图像利用Otsu算法算出一个阈值，将其作为一个高阈值TH，高阈值的一半作为低阈值TL。

2 算法步骤小结

使用一个高斯滤波器平滑输入图像。
计算梯度幅值图像和角度图像。
对梯度幅值图像进行非极大值抑制。
将非极大值抑制获得的图像利用Otsu算法确定双阈值。
使用双阈值处理和连通域分析来检测与连接边缘。

具体内容可参照冈萨雷斯《数字图像处理》

具体代码如下：

//确定一个点的坐标是否在图像内
bool checkInRang(int r, int c, int rows, int cols) {if (r >= 0 && r < rows && c >= 0 && c < cols)return true;elsereturn false;
}//从确定边缘点出发，延长边缘
void EdgePoint_Trace(cv::Mat& edgeMag_noMaxsup, cv::Mat& edge, unsigned TL, 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){float mag = edgeMag_noMaxsup.at<float>(r + i, c + j);if (checkInRang(r + i, c + j, rows, cols) && mag >= TL)EdgePoint_Trace(edgeMag_noMaxsup, edge, TL, r + i, c + j, rows, cols);}}}
}/********************************mian函数入口***************************************/
int main()
{string path = "F:\\NoteImage\\lena.jpg";Mat SrcImage = imread(path);if (!SrcImage.data) {std::cout << "Could not open or find the image" << std::endl;return -1;}cv::Mat grayImage, cannyImage;cvtColor(SrcImage, grayImage, COLOR_BGR2GRAY);//使图像连续并可导GaussianBlur(grayImage, grayImage, Size(3, 3), 0, 0);cv::Mat gx, gy;cv::Mat mag, angle;Sobel(grayImage, gx, CV_32F, 1, 0, 3);Sobel(grayImage, gy, CV_32F, 0, 1, 3);//计算梯度幅值和梯度的方向（角度）cv::cartToPolar(gx, gy, mag, angle, true);//定义全黑非极大值抑制图像cv::Mat Non_maxImage = cv::Mat::zeros(grayImage.size(), CV_32FC1);int height = grayImage.rows;int width = grayImage.cols;//获得非极大值抑制图像for (int i = 1; i < height - 1; ++i){for (int j = 1; j < width - 1; ++j){float g_angle = angle.at<float>(i, j);float K_mag = mag.at<float>(i, j);//梯度方向在垂直方向if ((g_angle <= 112.5 && g_angle > 67.5) || (g_angle <= 292.5 && g_angle > 247.5)){if (K_mag >= mag.at<float>(i - 1, j) && K_mag >= mag.at<float>(i + 1, j))Non_maxImage.at<float>(i, j) = K_mag;}//梯度方向在水平方向else if (g_angle <= 22.5 || g_angle > 337.5 || (g_angle <= 202.5 && g_angle > 157.5)){if (K_mag >= mag.at<float>(i, j - 1) && K_mag >= mag.at<float>(i, j + 1))Non_maxImage.at<float>(i, j) = K_mag;}//梯度方向在+45方向else if ((g_angle <= 67.5 && g_angle > 22.5) || (g_angle <= 247.5 && g_angle > 202.5)){if (K_mag >= mag.at<float>(i - 1, j - 1) && K_mag >= mag.at<float>(i + 1, j + 1))      Non_maxImage.at<float>(i, j) = K_mag;}//梯度方向在-45方向else if ((g_angle <= 337.5 && g_angle > 292.5) || (g_angle <= 157.5 && g_angle > 112.5)){if (K_mag >= mag.at<float>(i + 1, j - 1) && K_mag >= mag.at<float>(i - 1, j + 1))     Non_maxImage.at<float>(i, j) = K_mag;}}}//双阈值处理--根据Otsu算出的阈值确定为高阈值，取高阈值的一半记为低阈值unsigned TH = Otsu_threshold(Non_maxImage);unsigned TL = TH * 0.5;cv::Mat My_cannyImage = cv::Mat::zeros(grayImage.size(), grayImage.type());for (int i = 1; i < height - 1; ++i){for (int j = 1; j < width - 1; ++j){float K_mag = Non_maxImage.at<float>(i, j);//大于高阈值确定为边缘点if (K_mag > TH)EdgePoint_Trace(Non_maxImage, My_cannyImage, TL, i, j, height, width);else if (K_mag < TL)My_cannyImage.at<uchar>(i, j) = 0;}}//和OpenCV自带函数做对比Canny(grayImage, cannyImage, TH, TL, 3, true);imshow("src", My_cannyImage);cv::waitKey(0);return 0;

双阈值边缘连接处理要点采用了大佬的方法：canny算子边缘检测原理与实现

试验图例：

2022/5/11更新

以上采用递归的方式去实现，代码简洁，但是当图像太大，在某些编译环境下，会有栈溢出的风险，其次算法只需要判断当前梯度方向状态是水平、垂直还是对角，并不需要实际去计算实际的梯度角，基于此，通过阅读源码和查找资料，做出了一些改进。

算法思路在上文中已有简要说明，下面直接给出代码：

//5×5高斯滤波
cv::Mat _gaussian_filter(const cv::Mat& mat)
{cv::Mat matDouble;mat.convertTo(matDouble, CV_64FC1);cv::Mat kernel = (cv::Mat_<double>(5, 5) <<2, 4, 5, 4, 2,4, 9, 12, 9, 4,5, 12, 15, 12, 5,4, 9, 12, 9, 4,2, 4, 5, 4, 2);kernel = kernel / 159;cv::Mat resDouble;cv::filter2D(matDouble, resDouble, -1, kernel, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);cv::Mat res;resDouble.convertTo(res, CV_8UC1);return res;
}//对滤波后的图利用sobel计算梯度，通过梯度角的tan值与tan22.5进行一些比较获取梯度角所属分区
//angle = 0-> horizontal, 1 -> vertical, 2 -> diagonal
void _sobel_gradient(const cv::Mat& mat, cv::Mat& dx, cv::Mat& dy, cv::Mat& magnitudes, cv::Mat& angles,int apertureSize, bool L2gradient)
{CV_Assert(apertureSize == 3 || apertureSize == 5);double scale = 1.0;cv::Sobel(mat, dx, CV_16S, 1, 0, apertureSize, scale, cv::BORDER_REPLICATE);cv::Sobel(mat, dy, CV_16S, 0, 1, apertureSize, scale, cv::BORDER_REPLICATE);const int TAN225 = 13573;           //tan22.5 * 2^15(2 << 15)angles = cv::Mat(mat.size(), CV_8UC1);  // 0-> horizontal, 1 -> vertical, 2 -> diagonalmagnitudes = cv::Mat::zeros(mat.rows + 2, mat.cols + 2, CV_32SC1);cv::Mat magROI = cv::Mat(magnitudes, cv::Rect(1, 1, mat.cols, mat.rows));for (int i = 0; i < mat.rows; i++){for (int j = 0; j < mat.cols; j++){short xs = dx.ptr<short>(i)[j];short ys = dy.ptr<short>(i)[j];int x = (int)std::abs(xs);int y = (int)std::abs(ys) << 15;if (L2gradient) {//magROI.ptr<int>(i)[j] = int(xs) * xs + int(ys) * ys;magROI.ptr<int>(i)[j] = (int)std::sqrt(xs * xs + ys * ys);}else {magROI.ptr<int>(i)[j] = std::abs(int(xs)) + std::abs(int(ys));}int tan225x = x * TAN225;if (y < tan225x) {  // horizontalangles.ptr<uchar>(i)[j] = 0;}else{int tan675x = tan225x + (x << 16);if (y > tan675x) {  // verticalangles.ptr<uchar>(i)[j] = 1;}else {  // diagonalangles.ptr<uchar>(i)[j] = 2;}}}}
}//根据angles将梯度图进行非极大值抑制得到NMSImage，对其利用OTSU算法计算阈值，
//计算得到的阈值为高阈值high，低阈值取0.5*high
void _calculate_hysteresis_threshold_value(const cv::Mat& dx, const cv::Mat& dy, const cv::Mat& magnitudes,const cv::Mat& angles, cv::Mat& NMSImage, int& low, int& high)
{NMSImage = cv::Mat::zeros(magnitudes.size(), magnitudes.type());      //CV_32SC1for (int i = 0; i < dx.rows; ++i){int r = i + 1;for (int j = 0; j < dx.cols; ++j){int c = j + 1;int m = magnitudes.ptr<int>(r)[c];uchar angle = angles.ptr<uchar>(i)[j];if (angle == 0)           //horizontal{   if (m > magnitudes.ptr<int>(r)[c - 1] && m >= magnitudes.ptr<int>(r)[c + 1])NMSImage.ptr<int>(r)[c] = m;}else if (angle == 1)  //vertical{if (m > magnitudes.ptr<int>(r - 1)[c] && m >= magnitudes.ptr<int>(r + 1)[c])NMSImage.ptr<int>(r)[c] = m;}else if (angle == 2)   //diagonal{short xs = dx.ptr<short>(i)[j];short ys = dy.ptr<short>(i)[j];if ((xs > 0 && ys > 0) || (xs < 0 && ys < 0)){   //45 degreeif (m > magnitudes.ptr<int>(r - 1)[c - 1] && m > magnitudes.ptr<int>(r + 1)[c + 1])NMSImage.ptr<int>(r)[c] = m;}else{ //135 degreeif (m > magnitudes.ptr<int>(r - 1)[c + 1] && m > magnitudes.ptr<int>(r + 1)[c - 1])NMSImage.ptr<int>(r)[c] = m;}}}}//利用Otsu对非极大值抑制图像进行处理，将计算得到的阈值作为高阈值high, 低阈值取高阈值的0.5倍cv::normalize(NMSImage, NMSImage, 0, 255, cv::NORM_MINMAX);NMSImage.convertTo(NMSImage, CV_8UC1);cv::Mat temp;high = (int)cv::threshold(NMSImage, temp, 0, 255, cv::THRESH_OTSU);low = (int)(0.5 * high);
}//对非极大值抑制后的图根据高低阈值进行标记，当当前像素小于low，则标记为1，当当前像素大于low且大于high，则标记为2
//当大于low小于high时标记为0，并将标记为2的像素坐标压入队列
void _non_maximum_suppression(const cv::Mat& NMSImage, cv::Mat& map, std::deque<int>& mapIndicesX,std::deque<int>& mapIndicesY, int low, int high)
{// 0 -> the pixel may be edge// 1 -> the pixel is not edge// 2 -> the pixel is edgemap = cv::Mat::ones(NMSImage.size(), CV_8UC1);for (int i = 0; i < NMSImage.rows; ++i){for (int j = 0; j < NMSImage.cols; ++j){int m = NMSImage.ptr<uchar>(i)[j];               //nms -> CV_8UC1if (m > low){if (m > high){map.ptr<uchar>(i)[j] = 2;mapIndicesX.push_back(j);mapIndicesY.push_back(i);}elsemap.ptr<uchar>(i)[j] = 0;}}}
}//双阈值滞后处理：根据队列中的像素坐标，进行8领域边缘点寻找，即在map中与2相连的0均认作为边缘点
void _hysteresis_thresholding(std::deque<int>& mapIndicesX, std::deque<int>& mapIndicesY, cv::Mat& map)
{while (!mapIndicesX.empty()){int r = mapIndicesY.back();int c = mapIndicesX.back();//获取到边缘点之后要将其弹出mapIndicesX.pop_back();mapIndicesY.pop_back();// top leftif (map.ptr<uchar>(r - 1)[c - 1] == 0) {mapIndicesX.push_back(c - 1);mapIndicesY.push_back(r - 1);map.ptr<uchar>(r - 1)[c - 1] = 2;}// topif (map.ptr<uchar>(r - 1)[c] == 0) {mapIndicesX.push_back(c);mapIndicesY.push_back(r - 1);map.ptr<uchar>(r - 1)[c] = 2;}// top rightif (map.ptr<uchar>(r - 1)[c + 1] == 0) {mapIndicesX.push_back(c + 1);mapIndicesY.push_back(r - 1);map.ptr<uchar>(r - 1)[c + 1] = 2;}// leftif (map.ptr<uchar>(r)[c - 1] == 0) {mapIndicesX.push_back(c - 1);mapIndicesY.push_back(r);map.ptr<uchar>(r)[c - 1] = 2;}// rightif (map.ptr<uchar>(r)[c + 1] == 0) {mapIndicesX.push_back(c + 1);mapIndicesY.push_back(r);map.ptr<uchar>(r)[c + 1] = 2;}// bottom leftif (map.ptr<uchar>(r + 1)[c - 1] == 0) {mapIndicesX.push_back(c - 1);mapIndicesY.push_back(r + 1);map.ptr<uchar>(r + 1)[c - 1] = 2;}// bottomif (map.ptr<uchar>(r + 1)[c] == 0) {mapIndicesX.push_back(c);mapIndicesY.push_back(r + 1);map.ptr<uchar>(r + 1)[c] = 2;}// bottom rightif (map.ptr<uchar>(r + 1)[c + 1] == 0) {mapIndicesX.push_back(c + 1);mapIndicesY.push_back(r + 1);map.ptr<uchar>(r + 1)[c + 1] = 2;}}
}cv::Mat _get_canny_result(const cv::Mat& map)
{cv::Mat dst(map.rows - 2, map.cols - 2, CV_8UC1);for (int i = 0; i < dst.rows; i++) {for (int j = 0; j < dst.cols; j++) {dst.ptr<uchar>(i)[j] = (map.ptr<uchar>(i + 1)[j + 1] == 2 ? 255 : 0);}}return dst;
}/*--------函数封装---------*/
//自适应阈值canny plus版本
cv::Mat Adaptive_Canny(const cv::Mat& src, int apertureSize, bool L2gradient)
{CV_Assert(src.type() == CV_8UC1);CV_Assert(apertureSize == 3 || apertureSize == 5);cv::Mat gaussianSrc = _gaussian_filter(src);cv::Mat dx, dy, magnitudes, angles;_sobel_gradient(gaussianSrc, dx, dy, magnitudes, angles, apertureSize, L2gradient);//非极大值抑制计算高低阈值int low, high;cv::Mat NMSImage;_calculate_hysteresis_threshold_value(dx, dy, magnitudes, angles, NMSImage, low, high);cv::Mat map;std::deque<int> mapIndicesX, mapIndicesY;_non_maximum_suppression(NMSImage, map, mapIndicesX, mapIndicesY, low, high);_hysteresis_thresholding(mapIndicesX, mapIndicesY, map);cv::Mat dst = _get_canny_result(map);return dst;
}//*****************************测试代码******************************//int main()
{std::string path = "F:\\NoteImage\\手掌2.2.jpg";cv::Mat src = cv::imread(path, cv::IMREAD_GRAYSCALE);if (!src.data) {std::cout << "Could not open or find the image" << std::endl;return -1;}cv::Mat gaussianSrc = _gaussian_filter(src);int apertureSize = 3;bool L2gradient = true;cv::Mat dx, dy, magnitudes, angles;_sobel_gradient(gaussianSrc, dx, dy, magnitudes, angles, apertureSize, L2gradient);//非极大值抑制计算高低阈值int low, high;cv::Mat NMSImage;_calculate_hysteresis_threshold_value(dx, dy, magnitudes, angles, NMSImage, low, high);cv::Mat map;std::deque<int> mapIndicesX, mapIndicesY;_non_maximum_suppression(NMSImage, map, mapIndicesX, mapIndicesY, low, high);_hysteresis_thresholding(mapIndicesX, mapIndicesY, map);cv::Mat dst = _get_canny_result(map);//利用计算出来的low和high传入opencv Canny进行对比cv::Mat opencvCanny;cv::Canny(gaussianSrc, opencvCanny, low, high, apertureSize, L2gradient);cv::imshow("dst", dst);cv::imshow("opencvCanny", opencvCanny);cv::waitKey(0);return 0;
}

结果分析：

 原图

试验结果

很神奇的是，将计算得到的阈值传入cv::canny(),与本文算法获得的结果图比较，有很大的差异，至于为什么会出现这个差异，在算法逻辑上想了很久都没有找到问题，但这个差异对我来说有好的一面：对于示例图，手掌是我们想要的前景图，在边缘分析中应该尽可能的去除背景成分。很显然，基于以上假设，本文算法达到了更好的效果，不仅提取到了完整的前景细腻边缘，还去除了一些背景边缘，误打误撞吧，这在某些场景下，还是挺有用的。

参考资料：基于改进Canny算子的锂电池极片表面缺陷检测
参考代码： B站大佬up主