Canny算子是边缘检测算子中最常用的一种，是公认性能优良的一种算子，常被其它边缘检测算子作为标准算子进行优劣分析。
1. 原理
Canny算法基本可以分为3个步骤：平滑、梯度计算、基于梯度值及梯度方向的候选点过滤
（1）平滑
图像梯度的计算对噪声很敏感，因此必须首先对其进行低通滤波。在这里使用5*5的高斯滤波器，为了提高滤波效率，使用近似的高斯模板（图1-1）
图1-1. 5*5的高斯模板
（2）梯度计算
这一步主要通过计算每个像素x、y方向的梯度，获得该点的梯度值和梯度方向。在这里使用一阶边缘检测微分算子Sobel（图1-2）
图1-2. sobel边缘检测算子
对每个点，由公式
可得到其梯度值，由公式
可得到对应的梯度方向
可得到其梯度值，由公式可得到对应的梯度方向
具体处理的时候，将角度归并到0、45、90、135这4个方向
（3）候选点过滤
这一步主要包含2块内容：基于梯度方向的非最大值抑制、基于双阈值的候选点剔除及断点连接
基于梯度方向的非最大值抑制
非最大值抑制主要考虑这种情况：处于边缘上的点，其邻域内梯度方向上（正方向和负方向）的点的梯度响应应该比该点小
然而实际中由梯度方向得到的邻域点的坐标往往是浮点值，因此需要对该值进行处理，映射到整数值
在这里介绍3种处理办法：
1. 将角度近似处理到0,45,90,135这4个方向，这样只需要对当前点向左、右、上、下移动就可以了。这种方式计算简单，但误差较大
2. 对浮点坐标进行插值（图1-3）
图1-3. 通过g1,g2 和g3,g4的加权和分别得到两个浮点坐标的梯度值（图片来自http://blog.youkuaiyun.com/likezhaobin/article/details/6892176）
考虑图示的情况，有dTemp1 = g1*weight + g2*(1-weight)，dTemp2 = g3*weight + g4*(1-weight)。其中weight由x和y方向的梯度值确定：
3. 考虑浮点坐标周围的4个整形坐标，取其均值作为dTemp1和dTemp2（参考自http://blog.youkuaiyun.com/likezhaobin/article/details/6892629里mike07026 的评论），增加半径值可以过滤离散点，但因为忽略了旁边的点，边缘会有膨胀效果
基于双阈值的候选点剔除及断点连接
为了进一步过滤噪声的影响，非极大值抑制后，Canny使用一高一低两个阈值处理梯度响应矩阵。
1.对响应值大于高阈值的点，直接标记为边缘点
2.对响应值小于高阈值，但大于低阈值的点，考察其8邻域，如果存在已被标记的边缘点，则也标记该点为边缘点
关于高低阈值的设置，文章http://blog.youkuaiyun.com/likezhaobin/article/details/6892629统计边缘点的数目，将边缘点按梯度响应从小到大排列，取位于第79%的那个点对应的梯度响应值作为高阈值，取其一半为低阈值
2. 代码
typedef struct mySize
{
int height;
int width;
}_size;
#define _MAX(a,b) ((a)>(b)?(a):(b))
#define _MIN(a,b) ((a)<(b)?(a):(b))
void _edgeCanny(unsigned char* img, unsigned char* out, _size imgSize, int lowThreshold, int highThreshold);
template<typename T1>
void convolution(unsigned char* srcImg, T1* dstImg, _size imgSize, int* kernel, _size kerSize, bool div);
void calcGradAngle(int* mx, int* my, float* grad, float* angle, _size imgsize);
void depress(int* mx, int* my, float* grad, float* angle, unsigned char* out, _size imgSize);
void depress_(float* grad, float* angle, unsigned char* out, _size imgSize, int radius);
void doubleThreshold(unsigned char* src, float* grad, _size imgSize, int low, int high);
void getThresholds(unsigned char* src, float* grad, _size imgSize, int* high, int* low);
void connectLow(unsigned char* src, float* grad, _size imgSize, int y, int x, int low);
void main()
{
cv::Mat img = cv::imread("../file/einstein.jpg", 0);
unsigned char* imgBuffer = new unsigned char[img.rows*img.cols];
for(int i=0; i<img.rows; i++)
{
uchar* ptr = img.ptr<uchar>(i);
for(int j=0; j<img.cols; j++)
imgBuffer[i*img.cols+j] = ptr[j];
}
_size srcSize;
srcSize.height = img.rows;
srcSize.width = img.cols;
unsigned char* outBuffer = new unsigned char[img.rows*img.cols];
_edgeCanny(imgBuffer, outBuffer, srcSize, 0, 0);
cv::Mat colorImg(img.size(), CV_8UC1, cv::Scalar::all(0));
for(int i=0; i<img.rows; i++)
{
uchar* ptr = colorImg.ptr<uchar>(i);
for(int j=0; j<img.cols; j++)
ptr[j] = outBuffer[i*img.cols+j];
}
cv::namedWindow("show");
cv::imshow("show", colorImg);
cv::waitKey(0);
}
void _edgeCanny(unsigned char* img, unsigned char* out, _size imgSize, int lowThreshold, int highThreshold)
{
unsigned char* smoothImg = new unsigned char[imgSize.height*imgSize.width];
int smKernel[25] = {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};
_size smSize;
smSize.height = 2; smSize.width = 2;
convolution(img, smoothImg, imgSize, smKernel, smSize, true);
int XKernel[9] = {-1,0,1,-2,0,2,-1,0,1};
int YKernel[9] = {-1,-2,-1,0,0,0,1,2,1};
smSize.height = 1; smSize.width = 1;
int* Mx = new int[imgSize.height*imgSize.width];
int* My = new int[imgSize.height*imgSize.width];
convolution(smoothImg, Mx, imgSize, XKernel, smSize, false);
convolution(smoothImg, My, imgSize, YKernel, smSize, false);
delete[] smoothImg;
float* gradient = new float[imgSize.height*imgSize.width];
float* angle = new float[imgSize.height*imgSize.width];
calcGradAngle(Mx, My, gradient, angle, imgSize);
depress(Mx, My, gradient, angle, out, imgSize);
// depress_(gradient, angle, out, imgSize, 1);
delete[] Mx;
delete[] My;
doubleThreshold(out, gradient, imgSize, lowThreshold, highThreshold);
delete[] gradient;
delete[] angle;
}
template<typename T1>
void convolution(unsigned char* srcImg, T1* dstImg, _size imgSize, int* kernel, _size kerSize, bool div)
{
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
int sumValue = 0;
int sumKernel = 0;
int count = 0;
for(int m=-kerSize.height; m<=kerSize.height; m++)
{
for(int n=-kerSize.width; n<=kerSize.width; n++)
{
int indexX = j+n;
int indexY = i+m;
if(indexX>=0 && indexX<imgSize.width && indexY>=0 && indexY<imgSize.height)
{
sumValue += int(srcImg[indexY*imgSize.width+indexX])*kernel[count];
sumKernel += kernel[count];
}
count++;
}
}
if(div)
{
sumValue /= sumKernel;
dstImg[i*imgSize.width+j] = _MAX(_MIN(sumValue,255), 0);
}
else
dstImg[i*imgSize.width+j] = sumValue;
}
}
}
void calcGradAngle(int* mx, int* my, float* grad, float* angle, _size imgsize)
{
for(int i=0; i<imgsize.height; i++)
{
for(int j=0; j<imgsize.width; j++)
{
grad[i*imgsize.width+j] = sqrt(pow(float(mx[i*imgsize.width+j]),2) +
pow(float(my[i*imgsize.width+j]),2));
angle[i*imgsize.width+j] = atan2(float(my[i*imgsize.width+j]), float(mx[i*imgsize.width+j])) *
180/PI;
if(angle[i*imgsize.width+j]<0)
angle[i*imgsize.width+j] += 360;
}
}
}
void depress(int* mx, int* my, float* grad, float* angle, unsigned char* out, _size imgSize)
{
float g1, g2, g3, g4;
float weight, v12, v34;
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
int index = i*imgSize.width+j;
if(i==0 || i==imgSize.height-1 || j==0 || j==imgSize.width-1
|| grad[index] <= 0.001)
out[index] = 0;
else
{
if((angle[index]>=90 && angle[index]<135) ||
(angle[index]>=270 && angle[index]<315))
{
/ g1 g2 /
/ C /
/ g3 g4 /
g1 = grad[index-imgSize.width-1];
g2 = grad[index-imgSize.width];
g3 = grad[index+imgSize.width];
g4 = grad[index+imgSize.width+1];
weight = abs(mx[index]/float(my[index]));
v12 = g1*weight + g2*(1-weight);
v34 = g4*weight + g3*(1-weight);
}
else if((angle[index]>=135 && angle[index]<180) ||
(angle[index]>=315 && angle[index]<360))
{
/ g1 /
/ g2 C g3 /
/ g4 /
g1 = grad[index-imgSize.width-1];
g2 = grad[index-1];
g3 = grad[index+1];
g4 = grad[index-imgSize.width+1];
weight = abs(my[index]/float(mx[index]));
v12 = g1*weight + g2*(1-weight);
v34 = g4*weight + g3*(1-weight);
}
else if((angle[index]>=45 && angle[index]<90) ||
(angle[index]>=225 && angle[index]<270))
{
/ g1 g2 /
/ C /
/ g4 g3 /
g1 = grad[index-imgSize.width];
g2 = grad[index-imgSize.width+1];
g3 = grad[index+imgSize.width];
g4 = grad[index+imgSize.width-1];
weight = abs(mx[index]/float(my[index]));
v12 = g2*weight + g1*(1-weight);
v34 = g4*weight + g3*(1-weight);
}
else if((angle[index]>=0 && angle[index]<45) ||
(angle[index]>=180 && angle[index]<225))
{
/ g1 /
/ g4 C g2 /
/ g3 /
g1 = grad[index-imgSize.width+1];
g2 = grad[index+1];
g3 = grad[index+imgSize.width-1];
g4 = grad[index-1];
weight = abs(my[index]/float(mx[index]));
v12 = g2*weight + g1*(1-weight);
v34 = g4*weight + g3*(1-weight);
}
if(grad[index]>=v12 && grad[index]>=v34)
out[index] = 128;
else
out[index] = 0;
}
}
}
}
void depress_(float* grad, float* angle, unsigned char* out, _size imgSize, int radius)
{
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
int index = i*imgSize.width+j;
if(i==0 || i==imgSize.height-1 || j==0 || j==imgSize.width-1
|| grad[index] <= 0.001)
out[index] = 0;
else
{
float realX1 = j + radius*cos(angle[i*imgSize.width+j]);
float realY1 = i + radius*sin(angle[i*imgSize.width+j]);
float realX2 = j - radius*cos(angle[i*imgSize.width+j]);
float realY2 = i - radius*sin(angle[i*imgSize.width+j]);
int count = 0;
float sumValue1 = 0, sumValue2 = 0;
for(int m=int(realY1); m<=int(realY1+1); m++)
{
for(int n=int(realX1); n<=int(realX1+1); n++)
{
if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
{
sumValue1 += grad[m*imgSize.width+n];
count++;
}
}
}
sumValue1 /= count;
count = 0;
for(int m=int(realY2); m<=int(realY2+1); m++)
{
for(int n=int(realX2); n<=int(realX2+1); n++)
{
if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
{
sumValue2 += grad[m*imgSize.width+n];
count++;
}
}
}
sumValue2 /= count;
if(grad[index]>=sumValue1 && grad[index]>=sumValue2)
out[index] = 128;
else
out[index] = 0;
}
}
}
}
void doubleThreshold(unsigned char* src, float* grad, _size imgSize, int low, int high)
{
if(low<0 || high<0 || low>=high)
getThresholds(src, grad, imgSize, &high, &low);
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
if(int(src[i*imgSize.width+j]) == 128 && grad[i*imgSize.width+j] >= high)
{
src[i*imgSize.width+j] = 255;
connectLow(src, grad, imgSize, i, j, low);
}
}
}
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
if(int(src[i*imgSize.width+j]) != 255)
src[i*imgSize.width+j] = 0;
}
}
}
void getThresholds(unsigned char* src, float* grad, _size imgSize, int* high, int* low)
{
int hist[1200];
int Ehist[1200];
for(int i=0; i<1200; i++)
{
hist[i] = 0;
Ehist[i] = 0;
}
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
if(int(src[i*imgSize.width+j]) == 128)
hist[int(grad[i*imgSize.width+j]+0.5)]++;
}
}
for(int i=0; i<1200; i++)
Ehist[i] = Ehist[(i-1)<0?0:(i-1)] + hist[i];
int pos = int(Ehist[1199]*0.79+0.5);
for(int i=0; i<1200; i++)
{
if(Ehist[i]>=pos)
{
*high = i;
break;
}
}
*low = int(0.5*(*high)+0.5);
}
void connectLow(unsigned char* src, float* grad, _size imgSize, int y, int x, int low)
{
int xNum[8] = {1,1,0,-1,-1,-1,0,1};
int yNum[8] = {0,1,1,1,0,-1,-1,-1};
for(int i=0; i<8; i++)
{
int xx = x+xNum[i];
int yy = y+yNum[i];
if(xx>=0 && xx<imgSize.width && yy>=0 && yy<imgSize.height)
{
if(int(src[yy*imgSize.width+xx]) == 128 && grad[yy*imgSize.width+xx] >= low)
{
src[yy*imgSize.width+xx] = 255;
connectLow(src, grad, imgSize, yy, xx, low);
}
}
}
}效果（(a)为opencv函数结果(阈值取140,70), (b)为对浮点坐标插值的效果，(c)和(d)为在浮点坐标周围取4个整数值的效果，前者半径为1，后者为4）
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
