大家好，又见面了，我是你们的朋友全栈君。
文章目录
- Canny原理
- 非极大抑制
- 图像坐标系
- 边缘方向区间
- Canny算法的实现(1.0版)
- Canny算法的实现(2.0版)
- 实验结果
Canny原理
Canny的原理就不细说了，冈萨雷斯的《数字图像处理》（中文第三版）P463~465讲解的比较清楚，主要就四个步骤：
- 对图像进行高斯滤波(由于实际实现的时候使用Sobel计算梯度，Sobel具有滤波作用，所以实际的实现省略了高斯滤波)
- 计算梯度大小和梯度方向
- 对梯度幅值图像进行非极大抑制
- 双阈值处理和连接性分析（通常这一步与非极大抑制并行，详见下面的代码）
下面重点说一下非极大抑制。
非极大抑制
对一幅图像计算梯度大小和梯度方向后，需要进行非极大抑制，一般都是通过计算梯度方向，沿着梯度方向，判断该像素点的梯度大小是否是极大值。这里主要说一下方向的判断。
图像坐标系
OpenCV中采用下面的坐标系
《数字图像处理》（中文第三版）这本书中图像坐标顺序与上述坐标系是相反的，看书的时候，注意一下坐标的顺序。
我们这里采用的图像坐标系与OpenCV一致。
本文中采用顺时针角度为正（OpenCV中规定逆时针为正）
边缘方向区间
非极大抑制中，通常将边缘量化为4个方向，水平，垂直，45°和-45°，实际中，通过定义一个方向角方位，在该方位内认为是某一方向的边缘，实现中，我们通过计算梯度方向的范围从而判断边缘的方向（边缘的方向与梯度方向垂直）。
由于边缘方向没有正负，梯度方向在-22.5°到 22.5和157.5°到-157.5°表示的是同一个边缘，所以这里梯度方向 θ theta θ只考虑180°区间内。如上图所示：
- 当 − 22.5 ° < θ < 22.5 ° -22.5°<theta<22.5° −22.5°<θ<22.5°的时候，梯度方向为0°，也就是垂直边缘
- 当 22.5 ° < θ < 67.5 ° 22.5°<theta<67.5° 22.5°<θ<67.5°的时候，梯度方向为45°，也就是135°边缘
- 其他两个方向依次类推
实现中，通过三角函数的性质计算，比如第一种情况，使用下面条件来判断 − 22.5 ° < θ < 22.5 ° -22.5°<theta<22.5° −22.5°<θ<22.5° => t a n ( − 22.5 ° ) < t a n θ < t a n ( 22.5 ° ) tan(-22.5°)<tantheta<tan(22.5°) tan(−22.5°)<tanθ<tan(22.5°) => t a n ( − 22.5 ° ) < f y f x < t a n ( 22.5 ° ) tan(-22.5°)<{fy over fx } <tan(22.5°) tan(−22.5°)<fxfy<tan(22.5°) => f x ∗ t a n ( − 22.5 ° ) < f y < f x ∗ t a n ( 22.5 ° ) fx*tan(-22.5°)<{fy } <fx* tan(22.5°) fx∗tan(−22.5°)<fy<fx∗tan(22.5°) 其中 f x 和 f y fx和fy fx和fy表示x方向和y方向的偏导数。
其他三种情况读者可以自行推导。
判断出梯度方向后，就可以进行非极大值抑制了。还是以第一种情况为例，比如我计算出了P5这个像素点处的梯度方向为0°(180°)，则这个时候，我们要判断的条件就是if(M(P5)>=M(P4)&&M(P5)>=M(P6)),也就是P5处的梯度值是否是极大值，其中M（）表示该像素点的梯度大小。
Canny算法的实现(1.0版)
代码语言：javascript复制#define CANNY_SHIFT 16
#define TAN_225 (int)(0.4142135623730950488016887242097*(1 << CANNY_SHIFT));
#define TAN_675 (int)(2.4142135623730950488016887242097*(1 << CANNY_SHIFT));
void Canny1(const Mat &srcImage, Mat &dstImage, double lowThreshold, double highThreshold, int sizeOfAperture, bool L2)
{
// 只支持灰度图
CV_Assert(srcImage.type() == CV_8UC1);
dstImage.create(srcImage.size(), srcImage.type());
// L2范数计算边缘强度的时候，距离采用平方的方式，所以阈值也需要采用平方
if (L2)
{
lowThreshold = std::min(32767.0, lowThreshold);
highThreshold = std::min(32767.0, highThreshold);
if (lowThreshold > 0) lowThreshold *= lowThreshold;
if (highThreshold > 0) highThreshold *= highThreshold;
}
// 计算fx,fy,强度图
Mat fx(srcImage.size(), CV_32SC1);
Mat fy(srcImage.size(), CV_32SC1);
Mat enlargedImage;
Mat magnitudeImage(srcImage.rows 2, srcImage.cols 2, CV_32SC1);
magnitudeImage.setTo(Scalar(0));
copyMakeBorder(srcImage, enlargedImage, 1, 1, 1, 1, cv::BORDER_REPLICATE);
int stepOfEnlargedImage = enlargedImage.cols;
int stepOffx = fx.cols;
int height = srcImage.rows;
int width = srcImage.cols;
uchar *rowOfEnlargedImage = enlargedImage.data stepOfEnlargedImage 1;
int *rowOffx = (int *)fx.data;
int *rowOffy = (int *)fy.data;
int *rowOfMagnitudeImage = (int *)magnitudeImage.data stepOfEnlargedImage 1;
for (int y = 0; y <= height - 1; y, rowOfEnlargedImage = stepOfEnlargedImage, rowOfMagnitudeImage = stepOfEnlargedImage, rowOffx = stepOffx, rowOffy = stepOffx)
{
uchar *colOfEnlargedImage = rowOfEnlargedImage;
int *colOffx = rowOffx;
int *colOffy = rowOffy;
int *colOfMagnitudeImage = rowOfMagnitudeImage;
for (int x = 0; x <= width - 1; x, colOfEnlargedImage, colOffx, colOffy, colOfMagnitudeImage)
{
// fx
colOffx[0] = colOfEnlargedImage[1 - stepOfEnlargedImage] 2 * colOfEnlargedImage[1] colOfEnlargedImage[1 stepOfEnlargedImage] -
colOfEnlargedImage[-1 - stepOfEnlargedImage] - 2 * colOfEnlargedImage[-1] - colOfEnlargedImage[-1 stepOfEnlargedImage];
// fy
colOffy[0] = colOfEnlargedImage[stepOfEnlargedImage - 1] 2 * colOfEnlargedImage[stepOfEnlargedImage] colOfEnlargedImage[stepOfEnlargedImage 1] -
colOfEnlargedImage[-stepOfEnlargedImage - 1] - 2 * colOfEnlargedImage[-stepOfEnlargedImage] - colOfEnlargedImage[-stepOfEnlargedImage 1];
// 计算边缘强度，由于只是用于比较，为了加快速度，只计算平方和
if (L2)
{
colOfMagnitudeImage[0] = colOffx[0] * colOffx[0] colOffy[0] * colOffy[0];
}
else
{
colOfMagnitudeImage[0] = std::abs(colOffx[0]) std::abs(colOffy[0]);
}
}
}
// 非极大抑制，同时标记图做标记,双阈值处理
// 0 - 可能是边缘
// 1 - 不是边缘
// 2 - 一定是边缘
Mat labelImage(srcImage.rows 2, srcImage.cols 2, CV_8UC1);
memset(labelImage.data, 1, labelImage.rows*labelImage.cols);
rowOffx = (int *)fx.data;
rowOffy = (int *)fy.data;
rowOfMagnitudeImage = (int *)magnitudeImage.data stepOfEnlargedImage 1;
uchar *rowOfLabelImage = labelImage.data stepOfEnlargedImage 1;
queue<uchar*> queueOfEdgePixel;
for (int y = 0; y <= height - 1; y, rowOfMagnitudeImage = stepOfEnlargedImage, rowOffx = stepOffx, rowOffy = stepOffx, rowOfLabelImage = stepOfEnlargedImage)
{
int *colOffx = rowOffx;
int *colOffy = rowOffy;
int *colOfMagnitudeImage = rowOfMagnitudeImage;
uchar *colOfLabelImage = rowOfLabelImage;
for (int x = 0; x <= width - 1; x, colOffx, colOffy, colOfMagnitudeImage, colOfLabelImage)
{
int fx = colOffx[0];
int fy = colOffy[0];
// 该像素点才有可能是边缘点
if (colOfMagnitudeImage[0] > lowThreshold)
{
// 非极大抑制
fy = fy*(1 << CANNY_SHIFT);
int tan225 = fx * TAN_225;
int tan675 = fx * TAN_675;
// 梯度方向0
if (fy>-1 * tan225 && fy < tan225)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[-1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[1])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
// 梯度方向45
if (fy > tan225&&fy < tan675)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[stepOfEnlargedImage 1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage - 1])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
// 梯度方向90
if (fy>tan675 || fy<-tan675)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[stepOfEnlargedImage] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
// 梯度方向135
if (fy >= -1 * tan675&&fy <= -1 * tan225)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >=colOfMagnitudeImage[stepOfEnlargedImage - 1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage 1])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
}
}
}
// 连接性分析,这里采用队列实现
while (!queueOfEdgePixel.empty())
{
uchar *m = queueOfEdgePixel.front();
queueOfEdgePixel.pop();
// 在8领域搜索
if (!m[-1])
{
m[-1] = 2;
queueOfEdgePixel.push(m - 1);
}
if (!m[1])
{
m[1] = 2;
queueOfEdgePixel.push(m 1);
}
if (!m[-stepOfEnlargedImage - 1])
{
m[-stepOfEnlargedImage - 1] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage - 1);
}
if (!m[-stepOfEnlargedImage])
{
m[-stepOfEnlargedImage] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage);
}
if (!m[-stepOfEnlargedImage 1])
{
m[-stepOfEnlargedImage 1] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage 1);
}
if (!m[stepOfEnlargedImage - 1])
{
m[stepOfEnlargedImage - 1] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage - 1);
}
if (!m[stepOfEnlargedImage])
{
m[stepOfEnlargedImage] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage);
}
if (!m[stepOfEnlargedImage 1])
{
m[stepOfEnlargedImage 1] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage 1);
}
}
// 最后生成边缘图
rowOfLabelImage = labelImage.data stepOfEnlargedImage 1;
uchar *rowOfDst = dstImage.data;
for (int y = 0; y <= height - 1; y, rowOfLabelImage = stepOfEnlargedImage, rowOfDst = stepOffx)
{
uchar *colOfLabelImage = rowOfLabelImage;
uchar *colOfDst = rowOfDst;
for (int x = 0; x <= width - 1; x, colOfDst, colOfLabelImage)
{
if (colOfLabelImage[0] == 2)
colOfDst[0] = 255;
else
{
colOfDst[0] = 0;
}
}
}
}Canny算法的实现(2.0版)
后来发现非极大抑制中梯度方向区间的计算可以使用绝对值求解，而且采用绝对值求解后，检测效果要好于1.0版本，修改后程序如下
代码语言：javascript复制#define CANNY_SHIFT 16
#define TAN_225 (int)(0.4142135623730950488016887242097*(1 << CANNY_SHIFT));
#define TAN_675 (int)(2.4142135623730950488016887242097*(1 << CANNY_SHIFT));
// 转为绝对值求解
void Canny2(const Mat &srcImage, Mat &dstImage, double lowThreshold, double highThreshold, int sizeOfAperture, bool L2)
{
// 只支持灰度图
CV_Assert(srcImage.type() == CV_8UC1);
dstImage.create(srcImage.size(), srcImage.type());
// L2范数计算边缘强度的时候，距离采用平方的方式，所以阈值也需要采用平方
if (L2)
{
lowThreshold = std::min(32767.0, lowThreshold);
highThreshold = std::min(32767.0, highThreshold);
if (lowThreshold > 0) lowThreshold *= lowThreshold;
if (highThreshold > 0) highThreshold *= highThreshold;
}
// 计算fx,fy,强度图
Mat fx(srcImage.size(), CV_32SC1);
Mat fy(srcImage.size(), CV_32SC1);
Mat enlargedImage;
Mat magnitudeImage(srcImage.rows 2, srcImage.cols 2, CV_32SC1);
magnitudeImage.setTo(Scalar(0));
copyMakeBorder(srcImage, enlargedImage, 1, 1, 1, 1, cv::BORDER_REPLICATE);
int stepOfEnlargedImage = enlargedImage.cols;
int stepOffx = fx.cols;
int height = srcImage.rows;
int width = srcImage.cols;
uchar *rowOfEnlargedImage = enlargedImage.data stepOfEnlargedImage 1;
int *rowOffx = (int *)fx.data;
int *rowOffy = (int *)fy.data;
int *rowOfMagnitudeImage = (int *)magnitudeImage.data stepOfEnlargedImage 1;
for (int y = 0; y <= height - 1; y, rowOfEnlargedImage = stepOfEnlargedImage, rowOfMagnitudeImage = stepOfEnlargedImage, rowOffx = stepOffx, rowOffy = stepOffx)
{
uchar *colOfEnlargedImage = rowOfEnlargedImage;
int *colOffx = rowOffx;
int *colOffy = rowOffy;
int *colOfMagnitudeImage = rowOfMagnitudeImage;
for (int x = 0; x <= width - 1; x, colOfEnlargedImage, colOffx, colOffy, colOfMagnitudeImage)
{
// fx
colOffx[0] = colOfEnlargedImage[1 - stepOfEnlargedImage] 2 * colOfEnlargedImage[1] colOfEnlargedImage[1 stepOfEnlargedImage] -
colOfEnlargedImage[-1 - stepOfEnlargedImage] - 2 * colOfEnlargedImage[-1] - colOfEnlargedImage[-1 stepOfEnlargedImage];
// fy
colOffy[0] = colOfEnlargedImage[stepOfEnlargedImage - 1] 2 * colOfEnlargedImage[stepOfEnlargedImage] colOfEnlargedImage[stepOfEnlargedImage 1] -
colOfEnlargedImage[-stepOfEnlargedImage - 1] - 2 * colOfEnlargedImage[-stepOfEnlargedImage] - colOfEnlargedImage[-stepOfEnlargedImage 1];
// 计算边缘强度，由于只是用于比较，为了加快速度，只计算平方和
if (L2)
{
colOfMagnitudeImage[0] = colOffx[0] * colOffx[0] colOffy[0] * colOffy[0];
}
else
{
colOfMagnitudeImage[0] = std::abs(colOffx[0]) std::abs(colOffy[0]);
}
}
}
// 非极大抑制，同时标记图做标记,双阈值处理
// 0 - 可能是边缘
// 1 - 不是边缘
// 2 - 一定是边缘
Mat labelImage(srcImage.rows 2, srcImage.cols 2, CV_8UC1);
memset(labelImage.data, 1, labelImage.rows*labelImage.cols);
rowOffx = (int *)fx.data;
rowOffy = (int *)fy.data;
rowOfMagnitudeImage = (int *)magnitudeImage.data stepOfEnlargedImage 1;
uchar *rowOfLabelImage = labelImage.data stepOfEnlargedImage 1;
queue<uchar*> queueOfEdgePixel;
for (int y = 0; y <= height - 1; y, rowOfMagnitudeImage = stepOfEnlargedImage, rowOffx = stepOffx, rowOffy = stepOffx, rowOfLabelImage = stepOfEnlargedImage)
{
int *colOffx = rowOffx;
int *colOffy = rowOffy;
int *colOfMagnitudeImage = rowOfMagnitudeImage;
uchar *colOfLabelImage = rowOfLabelImage;
for (int x = 0; x <= width - 1; x, colOffx, colOffy, colOfMagnitudeImage, colOfLabelImage)
{
int fx = colOffx[0];
int fy = colOffy[0];
int abs_fx = std::abs(fx);
int abs_fy = std::abs(fy);
// 该像素点才有可能是边缘点
if (colOfMagnitudeImage[0] > lowThreshold)
{
// 非极大抑制
abs_fy = abs_fy << CANNY_SHIFT;
int tan225 = abs_fx * TAN_225;
int tan675 = abs_fx * TAN_675;
// 梯度方向0°
if (abs_fy< tan225)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[-1] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[1])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
// 梯度方向90°
if (abs_fy>tan675)
{
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[stepOfEnlargedImage] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
// 梯度方向45°/-45°
if (abs_fy > tan225&&abs_fy<tan675)
{
int s = (fy*fx) < 0 ? -1 : 1;
// 极大值，有可能是边缘
if (colOfMagnitudeImage[0] >= colOfMagnitudeImage[-stepOfEnlargedImage - s] && colOfMagnitudeImage[0] >= colOfMagnitudeImage[stepOfEnlargedImage s])
{
// 大于高阈值，是边缘，标记为2
if (colOfMagnitudeImage[0] > highThreshold)
{
// 进入队列，并设置标记
colOfLabelImage[0] = 2;
queueOfEdgePixel.push(colOfLabelImage);
}
else
{
// 有可能是边缘，标记为0
colOfLabelImage[0] = 0;
}
}
}
}
}
}
// 连接性分析,这里采用队列实现
while (!queueOfEdgePixel.empty())
{
uchar *m = queueOfEdgePixel.front();
queueOfEdgePixel.pop();
// 在8领域搜索
if (!m[-1])
{
m[-1] = 2;
queueOfEdgePixel.push(m - 1);
}
if (!m[1])
{
m[1] = 2;
queueOfEdgePixel.push(m 1);
}
if (!m[-stepOfEnlargedImage - 1])
{
m[-stepOfEnlargedImage - 1] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage - 1);
}
if (!m[-stepOfEnlargedImage])
{
m[-stepOfEnlargedImage] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage);
}
if (!m[-stepOfEnlargedImage 1])
{
m[-stepOfEnlargedImage 1] = 2;
queueOfEdgePixel.push(m - stepOfEnlargedImage 1);
}
if (!m[stepOfEnlargedImage - 1])
{
m[stepOfEnlargedImage - 1] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage - 1);
}
if (!m[stepOfEnlargedImage])
{
m[stepOfEnlargedImage] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage);
}
if (!m[stepOfEnlargedImage 1])
{
m[stepOfEnlargedImage 1] = 2;
queueOfEdgePixel.push(m stepOfEnlargedImage 1);
}
}
// 最后生成边缘图
rowOfLabelImage = labelImage.data stepOfEnlargedImage 1;
uchar *rowOfDst = dstImage.data;
for (int y = 0; y <= height - 1; y, rowOfLabelImage = stepOfEnlargedImage, rowOfDst = stepOffx)
{
uchar *colOfLabelImage = rowOfLabelImage;
uchar *colOfDst = rowOfDst;
for (int x = 0; x <= width - 1; x, colOfDst, colOfLabelImage)
{
if (colOfLabelImage[0] == 2)
colOfDst[0] = 255;
else
{
colOfDst[0] = 0;
}
}
}
}实验结果
实验代码
代码语言：javascript复制int main(int argc, char *argv[])
{
Mat srcImage = imread("D:/Image/Gray/Lena512.bmp", -1);
Mat canny1,canny,canny2;
Canny1(srcImage, canny1, 50, 150, 3, false);
Canny2(srcImage, canny3, 50, 150, 3, false);
Canny(srcImage, canny, 50, 150);
imwrite("D:/Canny.bmp", canny);
imwrite("D:/Canny1.bmp", canny1);
imwrite("D:/Canny2.bmp", canny2);
return 0;
}使用的是标准的Lena图
OpenCV 的处理结果为：
Canny1 的结果:
Canny2 的结果:
结果可以看出，Canny1丢失了一些垂直边缘，改进后的Canny2与OpenCV处理结果基本一致。
完整工程见github项目：QQImageProcess_OpenCV 其中Canny边缘检测的实现在 Src/ImageProcess/Edge.h中的Canny_系列函数中。
