学习记录…
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;
else
return 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 -> diagonal
magnitudes = 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) { // horizontal
angles.ptr<uchar>(i)[j] = 0;
}
else
{
int tan675x = tan225x + (x << 16);
if (y > tan675x) { // vertical
angles.ptr<uchar>(i)[j] = 1;
}
else { // diagonal
angles.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_32SC1
for (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 degree
if (m > magnitudes.ptr<int>(r - 1)[c - 1] && m > magnitudes.ptr<int>(r + 1)[c + 1])
NMSImage.ptr<int>(r)[c] = m;
}
else
{ //135 degree
if (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 edge
map = 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_8UC1
if (m > low)
{
if (m > high)
{
map.ptr<uchar>(i)[j] = 2;
mapIndicesX.push_back(j);
mapIndicesY.push_back(i);
}
else
map.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 left
if (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;
}
// top
if (map.ptr<uchar>(r - 1)[c] == 0)
{
mapIndicesX.push_back(c);
mapIndicesY.push_back(r - 1);
map.ptr<uchar>(r - 1)[c] = 2;
}
// top right
if (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;
}
// left
if (map.ptr<uchar>(r)[c - 1] == 0)
{
mapIndicesX.push_back(c - 1);
mapIndicesY.push_back(r);
map.ptr<uchar>(r)[c - 1] = 2;
}
// right
if (map.ptr<uchar>(r)[c + 1] == 0)
{
mapIndicesX.push_back(c + 1);
mapIndicesY.push_back(r);
map.ptr<uchar>(r)[c + 1] = 2;
}
// bottom left
if (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;
}
// bottom
if (map.ptr<uchar>(r + 1)[c] == 0)
{
mapIndicesX.push_back(c);
mapIndicesY.push_back(r + 1);
map.ptr<uchar>(r + 1)[c] = 2;
}
// bottom right
if (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主
扫一扫
417
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
