网上有很多关于Canny实现的文章,但是多数在最后的边缘追踪时采用递归实现或者采用非递归实现个人认为效果又不好或者修改领域跟踪范围复杂。
采用递归实现可以参考:http://www.math.tau.ac.il/~turkel/notes/canny_source.c
matlab实现可以参考http://justin-liang.com/tutorials/canny/中给的源码,我就是根据此链接将matlab实现转换成C实现。
一、基本原理:
1、RGB转换成灰度:
上面是彩色图,下面是 灰度图。
2、高斯模糊。主要是滤噪,噪声可能影响生成边缘,所以这一步一般不可少。高斯模糊核的生成,首先根据sigma的值求得核大小,然后根据高斯公式计算卷积核。
kernel_size = 2 * ceil(2 * sigma) + 1;这里采用的是matlab中的imgaussfilt实现方式。冈萨雷斯的《数字图像处理》中建议kernel_size为3倍sigma的值。
高斯公式如下:
f
(
x
,
y
)
=
A
exp
(
−
(
x
−
x
o
)
2
+
(
y
−
y
o
)
2
2
σ
2
)
{\displaystyle f(x,y)=A\exp \left(-{\frac {(x-x_{o})^{2}+(y-y_{o})^{2}}{2\sigma ^{2}}}\right)}
f(x,y)=Aexp(−2σ2(x−xo)2+(y−yo)2)
原点为(0,0),在计算时常数A可以不用考虑,因为在生成高斯模糊核时还需要归一化处理。
3、使用sobel算子分别求取x方向和y方向梯度,x方向和y方向sobel算子如下:
G
x
=
[
+
1
0
−
1
+
2
0
−
2
+
1
0
−
1
]
G
y
=
[
+
1
+
2
+
1
0
0
0
−
1
−
2
−
1
]
{\displaystyle {G}_{x}={\begin{bmatrix}+1&0&-1\\+2&0&-2\\+1&0&-1\end{bmatrix}} \quad \quad{G} _{y}={\begin{bmatrix}+1&+2&+1\\0&0&0\\-1&-2&-1\end{bmatrix}} }
Gx=⎣⎡+1+2+1000−1−2−1⎦⎤Gy=⎣⎡+10−1+20−2+10−1⎦⎤
x方向的梯度图:
y方向的梯度图:
求取梯度幅值和方向,公式如下:
G
=
G
x
2
+
G
y
2
\mathbf {G} ={\sqrt {{\mathbf {G} _{x}}^{2}+{\mathbf {G} _{y}}^{2}}}
G=Gx2+Gy2
Θ
=
atan
(
G
y
G
x
)
{\displaystyle \mathbf {\Theta } =\operatorname {atan} \left({\mathbf {G} _{y} \over \mathbf {G} _{x}}\right)}
Θ=atan(GxGy)
幅值图:
4、非最大值抑制。由幅值图可以看出,梯度图主要是一些山脊线(ridges)组成,比较粗的,而边缘一般是比较细的,一个像素宽度大小,所以需要进行细化这些山脊线。非最大值抑制只是细化山脊线的一种方法,这种方法的本质在于指定梯度法向的离散数,也就是把360度方向离散化为几个方向。一般采用3x3领域,离散出4个方向,如下图:
我们知道,边缘方向与梯度方向垂直,也即是梯度方向与边缘法向平行,例如下图中18度方向的边缘:
对于3x3领域,设d1,d2,d3,d4分别表示水平、-45度、垂直、45度基础边缘方向,对角度图中的每一点a(x,y),找到最接近的基础边缘方向d,设K为(x,y)处的幅值,如果K小于d方向上任意一个幅值,则该处(x,y)幅值为0,否则不变。
非最大值移植结果图:
可以看出,边缘比较细了。
5、双阈值处理。从非最大值抑制的结果可以看出,有的边缘比较明显,如图中的比较亮的部分,而又的边缘比较暗,如图中相对暗的边缘,如果采用单阈值处理,则可能漏掉很多有用的边缘,如果采用3阈值以上处理,过程相对复杂,在复杂与效果的权衡下,采用双阈值处理,一个相对较高的阈值,得到较亮的边缘,这里成为第一类边缘点,一个较低阈值,滤掉不是不需要的边缘,剩下的可能是边缘也可能不是,判断的标准是如果这些可能的边缘点跟高阈值过滤得到的边缘点连接,则认为是边缘点,成为第二次类边缘点,否则不是。在第一次得到的边缘点与第二次类边缘点之间可能连接了多个第二次类边缘点,所以需要采用边缘追踪找出所有的与第一类边缘点连接的第二类边缘点,一般采用8领域的边缘点追踪算法,采用递归实现。
双阈值处理结果:
图中亮的边缘为第一类边缘点,灰的边缘可能是第二类边缘点,也可能不是。
采用边缘追踪并且去掉假边缘结果如下:
我采用的是5x5领域边缘追踪。
代码实现:
1、数据结构准备:
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <math.h>
typedef struct real_2d_array_t
{
float** arr;
float* data;
int rows;
int cols;
}real_2d;
typedef struct real_uint8_array_t
{
uint8_t** arr;
uint8_t* data;
int rows;
int cols;
}uint8_2d;
typedef struct int_vector_t
{
float* data;
int size;
}int_vec;
real_2d* create_real_2d(int rows, int cols);
void destroy_real_2d(real_2d** rd);
uint8_2d* create_uint8_2d(int rows, int cols);
void destroy_uint8_2d(uint8_2d** rd);
int_vec* create_int_vec(int size);
void destroy_int_vec(int_vec** rv);
real_2d* uint8_to_real(uint8_t* img, int w, int h);
uint8_2d* real_to_uint8(real_2d* img);
float get_real_2d_max(real_2d* rd);
void normalize_real_2d(real_2d* rd, float max_val);
void normalized_to_unnormalize(real_2d* rd);
这里定义两种二维数组,一种float型,一种uint8_t型,之所以定义2维数组,主要是方便访问图像像素值,把图像的尺寸也进行封装主要是让函数接口封装更加简洁。
相关实现:
#define FLT_MAX 3.402823e+38
#define FLT_MIN 1.175494e-38
real_2d* create_real_2d(int rows, int cols)
{
real_2d* rd = NULL;
rd = (real_2d*)calloc(1, sizeof(real_2d));
if (!rd)
{
return NULL;
}
rd->data = (float*)calloc(1, rows * cols * sizeof(float));
if (!rd->data)
{
destroy_real_2d(&rd);
return NULL;
}
rd->arr = (float**)calloc(1, rows * sizeof(float*));
if (!rd->arr)
{
destroy_real_2d(&rd);
return NULL;
}
int i;
for (i = 0; i < rows; i++)
{
rd->arr[i] = &rd->data[i*cols];
}
rd->cols = cols;
rd->rows = rows;
return rd;
}
void destroy_real_2d(real_2d** _rd)
{
if (!_rd) return;
real_2d* rd = *_rd;
if (!rd)
{
return;
}
if (rd->data)
{
free(rd->data);
rd->data = NULL;
}
if (rd->arr)
{
free(rd->arr);
rd->arr = NULL;
}
free(rd);
*_rd = NULL;
}
uint8_2d* create_uint8_2d(int rows, int cols)
{
uint8_2d* ud = NULL;
ud = (uint8_2d*)calloc(1, sizeof(uint8_2d));
if (!ud)
{
return NULL;
}
ud->data = (uint8_t*)calloc(1, rows * cols * sizeof(uint8_t));
if (!ud->data)
{
destroy_real_2d(&ud);
return NULL;
}
ud->arr = (uint8_t**)calloc(1, rows * sizeof(uint8_t*));
if (!ud->arr)
{
destroy_real_2d(&ud);
return NULL;
}
int i;
for (i = 0; i < rows; i++)
{
ud->arr[i] = &ud->data[i*cols];
}
ud->cols = cols;
ud->rows = rows;
return ud;
}
void destroy_uint8_2d(uint8_2d** _ud)
{
if (!_ud) return;
uint8_2d* ud = *_ud;
if (!ud)
{
return;
}
if (ud->data)
{
free(ud->data);
ud->data = NULL;
}
if (ud->arr)
{
free(ud->arr);
ud->arr = NULL;
}
free(ud);
*_ud = NULL;
}
int_vec* create_int_vec(int size)
{
if (size < 0)return NULL;
int_vec* iv = (int_vec*)calloc(1, sizeof(int_vec));
if (!iv) return iv;
iv->data = (int*)calloc(1, size*sizeof(int));
if (!iv->data)
{
destroy_int_vec(&iv);
return NULL;
}
iv->size = size;
return iv;
}
void destroy_int_vec(int_vec** _iv)
{
if (!_iv)return;
int_vec* iv = *_iv;
if (!iv) return;
if (iv->data)
{
free(iv->data);
iv->data = NULL;
}
free(iv);
*_iv = NULL;
}
real_2d* uint8_to_real(uint8_t* img, int w, int h)
{
real_2d* res = create_real_2d(h, w);
if (!res)return res;
int i, j;
float** arr = res->arr;
for (i = 0; i < h; i++)
{
for (j = 0; j < w; j++)
{
arr[i][j] = img[i*w + j];
}
}
res->cols = w;
res->rows = h;
return res;
}
uint8_2d* real_to_uint8(real_2d* img)
{
uint8_2d* res = NULL;
if (!img)return NULL;
int i, j;
int rows, cols;
int val;
rows = img->rows;
cols = img->cols;
res = create_uint8_2d(rows, cols);
if (!res)return NULL;
float** img_arr = img->arr;
uint8_t** res_arr = res->arr;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
val = round(img_arr[i][j]);
if (val > 255)val = 255;
if (val < 0) val = 0;
res_arr[i][j] = val;
}
}
return res;
}
float get_real_2d_max(real_2d* rd)
{
if (!rd)return FLT_MAX;
int i, j;
float max_val = FLT_MIN;
float** rd_arr = rd->arr;
int rows, cols;
rows = rd->rows;
cols = rd->cols;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
if (rd_arr[i][j] > max_val)
max_val = rd_arr[i][j];
}
}
return max_val;
}
void normalize_real_2d(real_2d* rd, float max_val)
{
if (!rd || max_val == 0.0)return;
int i, j;
float** rd_arr = rd->arr;
int rows, cols;
rows = rd->rows;
cols = rd->cols;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
rd_arr[i][j] /= max_val;
}
}
}
void normalized_to_unnormalize(real_2d* rd)
{
if (!rd)return;
int i, j;
float** rd_arr = rd->arr;
int rows, cols;
rows = rd->rows;
cols = rd->cols;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
rd_arr[i][j] *= 255.0;
}
}
}
2、canny相关方法:
高斯核生成:
real_2d* create_gaussian_kernel(float sigma)
{
int kernel_size;
kernel_size = 2 * ceil(2 * sigma) + 1;
real_2d* rd = create_real_2d(kernel_size, kernel_size);
if (!rd)return rd;
int i, j;
int half_size = kernel_size / 2;
float** rd_arr = rd->arr;
float sum = 0.0;
float sigma2 = sigma*sigma * 2;
float val = 0.0;
for (i = -half_size; i <= half_size; i++)
{
for (j = -half_size; j <= half_size; j++)
{
val = -1 * (i*i + j*j) / sigma2;
val = exp(val);
sum += val;
rd_arr[i + half_size][j + half_size] = val;
}
}
assert(sum != 0.0);
for (i = 0; i < kernel_size; i++)
{
for (j = 0; j < kernel_size; j++)
{
rd_arr[i][j] /= sum;
printf("%.4f ", rd_arr[i][j]);
}
printf("\n");
}
return rd;
}
高斯模糊:
void gaussian_blur(real_2d* img, real_2d* kernel, int r)
{
if (!img || !kernel)return;
int i, j;
int m, n;
int bound;
float sum = 0.0f;
int w = img->cols;
int h = img->rows;
real_2d* rd = create_real_2d(h, w);
if (!rd)return;
float** dst = rd->arr;
float** src = img->arr;
float** gc = kernel->arr;
int col, row;
bound = kernel->rows / 2;
for (i = 0; i < h; i++)
{
for (j = 0; j < w; j++)
{
sum = 0.0;
for (m = -bound; m <= bound; m++)
{
for (n = -bound; n <= bound; n++)
{
//if (i + m < 0)sum += 0;
//else if (i + m >= h)sum += 0;
//else if (j + n < 0)sum += 0;
//else if (j + n >= w)sum += 0;
//else sum += gc[m + bound][n + bound] * src[i+m][j+n];
if (i + m < 0)row=0;
else if (i + m >= h)row = h - 1;
else row = m + i;
if (j + n < 0)col = 0;
else if (j + n >= w)col = w - 1;
else col = j+n;
sum += gc[m + bound][n + bound] * src[row][col];
}
}
if (r)dst[i][j] = round(sum);
else dst[i][j] = sum;
}
}
for (i = 0; i < h; i++)
{
for (j = 0; j < w; j++)
{
src[i][j] = dst[i][j];
}
}
destroy_real_2d(&rd);
}
sobel算子求梯度:
real_2d* sobel_grad(real_2d* img, float sobel[SOL_COE_SIZE][SOL_COE_SIZE])
{
if (!img)return NULL;
int i, j;
int m, n;
int bound = SOL_COE_SIZE / 2;
float sum = 0.0f;
int w = img->cols;
int h = img->rows;
real_2d* rd = create_real_2d(h, w);
if (!rd)return NULL;
float** dst = rd->arr;
float** src = img->arr;
int row, col;
for (i = 0; i < h; i++)
{
for (j = 0; j < w; j++)
{
sum = 0.0;
for (m = -bound; m <= bound; m++)
{
for (n = -bound; n <= bound; n++)
{
if (i + m < 0)sum += 0;
else if (i + m >= h)sum += 0;
else if (j + n < 0)sum += 0;
else if (j + n >= w)sum += 0;
else sum += sobel[m + bound][n + bound] * src[i + m][j + n];
}
}
dst[i][j] = sum;
}
}
return rd;
}
计算幅值和角度,其实可以梯度放在一起:
void calculate_angle_and_mag(real_2d* sgx, real_2d* sgy, real_2d** dir, real_2d** mag)
{
if (!sgx || !sgy)return;
int rows;
int cols;
int i, j;
assert(sgx->rows == sgy->rows && sgx->cols == sgy->cols);
rows = sgx->rows;
cols = sgx->cols;
real_2d* direct = create_real_2d(rows, cols);
real_2d* magnitude = create_real_2d(rows, cols);
if (!direct || !magnitude)
{
destroy_real_2d(&direct);
destroy_real_2d(&magnitude);
return;
}
float** dir_arr = direct->arr;
float y;
float x;
float angle;
float** mag_arr = magnitude->arr;
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
y = sgy->arr[i][j];
x = sgx->arr[i][j];
angle = atan2(y, x) * 180.0 / M_PI;
dir_arr[i][j] = angle;
mag_arr[i][j] = sqrt(x*x + y*y);
}
}
*dir = direct;
*mag = magnitude;
}
非最大值抑制,同时求最大幅值,用作归一化:
void non_maximum_supression(real_2d* dir, real_2d* mag, real_2d** nms, float* mag_max)
{
if (!dir || !mag)return;
int rows;
int cols;
int i, j;
assert(dir->rows == mag->rows && dir->cols == mag->cols);
rows = mag->rows;
cols = mag->cols;
real_2d* res = create_real_2d(rows, cols);
if (!res)return NULL;
float** res_arr = res->arr;
float** dir_arr = dir->arr;
float** mag_arr = mag->arr;
float angle;
float mag_cur;
float mag1_max = -1.0;
int cur_max;
rows--;
cols--;
for (i = 1; i < rows; i++)
{
for (j = 1; j < cols; j++)
{
angle = dir_arr[i][j];
mag_cur = mag_arr[i][j];
cur_max = 0;
if ((angle > -22.5 && angle <= 22.5) || (angle < -157.5 && angle >= -180) || (angle > 157.5 && angle <= 180))
{
if (mag_cur >= mag_arr[i][j+1] && mag_cur >= mag_arr[i][j-1])
{
res_arr[i][j] = mag_arr[i][j];
}
else
{
res_arr[i][j] = 0;
}
}
else if ((angle > 22.5 && angle <= 67.5) || (angle < -112.5 && angle >= -157.5))
{
if (mag_cur >= mag_arr[i+1][j + 1] && mag_cur >= mag_arr[i-1][j - 1])
{
res_arr[i][j] = mag_arr[i][j];
}
else
{
res_arr[i][j] = 0;
}
}
else if ((angle > 67.5 && angle <= 112.5) || (angle >= -112.5 && angle <= -67.5))
{
if (mag_cur >= mag_arr[i + 1][j] && mag_cur >= mag_arr[i - 1][j])
{
res_arr[i][j] = mag_arr[i][j];
}
else
{
res_arr[i][j] = 0;
}
}
else if ((angle > 112.5 && angle <= 157.5) || (angle <= -22.5 && angle > -67.5))
{
if (mag_cur >= mag_arr[i + 1][j-1] && mag_cur >= mag_arr[i - 1][j+1])
{
res_arr[i][j] = mag_arr[i][j];
}
else
{
res_arr[i][j] = 0;
}
}
if (res_arr[i][j] > mag1_max)
mag1_max = res_arr[i][j];
}
}
*nms = res;
*mag_max = mag1_max;
}
双阈值处理,返回第一类边缘点和第二类边缘点的行号和列号:
void double_thresholding(real_2d* mag, float low_thres, float high_thres,
int_vec** _srv, int_vec** _scv, int* _scount,
int_vec** _wrv, int_vec** _wcv, int* _wcount)
{
if (!mag)return;
int rows, cols;
int i, j;
int scount;
int wcount;
float** mag_arr;
int_vec* srv;
int_vec* scv;
int_vec* wrv;
int_vec* wcv;
float mag_val;
rows = mag->rows;
cols = mag->cols;
mag_arr = mag->arr;
scount = 0;
wcount = 0;
srv = create_int_vec(rows*cols);
scv = create_int_vec(rows*cols);
wrv = create_int_vec(rows*cols);
wcv = create_int_vec(rows*cols);
if (!srv || !scv || !wrv || !wcv)
{
destroy_int_vec(&srv);
destroy_int_vec(&scv);
destroy_int_vec(&wrv);
destroy_int_vec(&wcv);
return;
}
for (i = 0; i < rows; i++)
{
for (j = 0; j < cols; j++)
{
mag_val = mag_arr[i][j];
if (mag_val > high_thres)
{
srv->data[scount] = i;
scv->data[scount] = j;
scount++;
mag_arr[i][j] = 1.0;
}
else if (mag_val < low_thres)
{
mag_arr[i][j] = 0;
}
else
{
wrv->data[wcount] = i;
wcv->data[wcount] = j;
wcount++;
}
}
}
*_scount = scount;
*_wcount = wcount;
*_scv = scv;
*_srv = srv;
*_wcv = wcv;
*_wrv = wrv;
}
边缘追踪,这里给出递归实现:
#define CHECK_BOUNDARY (row + i >= 0 && col + j >= 0 && row + i < rows && col + j < cols)
void find_connect_weak_edge(real_2d* mag, int row, int col)
{
if (!mag)return;
int rows, cols;
int i, j;
float val;
rows = mag->rows;
cols = mag->cols;
call_count++;
for (i = -2; i <= 2;i++)
{
for (j = -2; j <= 2; j++)
{
if CHECK_BOUNDARY
{
val = mag->arr[row + i][col + j];
if (val > 0.0 && val < 1)
{
mag->arr[row+i][col+j] = 1;
find_connect_weak_edge(mag, row + i, col + j);
}
}
}
}
}
void edge_track(real_2d* mag, int_vec* srv, int_vec* scv, int scount)
{
if (!mag || !srv || !scv)return;
int i;
for (i = 0; i < scount; i++)
{
find_connect_weak_edge(mag, srv->data[i], scv->data[i]);
}
}
非递归实现版本,这里的思想是,边缘点一定幅值大于较小阈值的点,但是大于较小阈值的点不一定是边缘点,需要与第一类边缘点直接在某邻域连接,这里使用5x5,所以,只要与第一类边缘点邻域连接点都是边缘点,也即是第二类点,然后把这些点加到第一类点的队列末尾,变成第一个类边缘点,同时第一类点个数增加1,当与第一类边缘点邻域相连结束,即完成边界追踪。这样实现需要改成7x7邻域,也很简单,只需要把下面循环中的2改为3即可。
void edge_track2(real_2d* mag, int_vec* srv, int_vec* scv, int scount)
{
if (!mag || !srv || !scv)return;
int i;
int m, n;
int rows, cols;
float val;
int row, col;
rows = mag->rows;
cols = mag->cols;
for (i = 0; i < scount; i++)
{
row = srv->data[i];
col = scv->data[i];
for (m = -2; m <= 2; m++)
{
for (n = -2; n <= 2; n++)
{
if (m + row >= 0 && n + col >= 0 && m + row < rows && n + col < cols)
{
val = mag->arr[row + m][col + n];
if (val > 0.0 && val < 1)
{
assert(scount < rows * cols);
mag->arr[row + m][col + n] = 1;
srv->data[scount] = row + m;
scv->data[scount] = col + n;
mag->arr[row+m][col+n] = 1;
scount++;
}
}
}
}
}
}
清理假边缘点:
void cleanup_weak_edge(real_2d* mag, int_vec* wrv, int_vec* wcv, int wcount)
{
if (!mag || !wrv || !wcv)return;
int i;
int row, col;
for (i = 0; i < wcount; i++)
{
row = wrv->data[i];
col = wcv->data[i];
if (mag->arr[row][col] != 1.0)
mag->arr[row][col] = 0.0;
}
}
canny边缘检测主函数:
void canny(uint8_t* data, int w, int h, uint8_t** edge)
{
real_2d* kernel = NULL;
real_2d* sgx = NULL;
real_2d* sgy = NULL;
real_2d* img = uint8_to_real(data, w, h);
real_2d* dir = NULL;
real_2d* mag = NULL;
real_2d* nms = NULL;
uint8_2d* e = NULL;
float sigma = 1.4;
float mag_max = 0.0;
float low_thred_ratio = 0.15;
float high_thred_ratio = 0.2;
float low_thred, high_thred;
float sobelx[3][3] = { { 1, 0, -1}, { 2, 0, -2}, { 1, 0, -1 } };
float sobely[3][3] = { { 1, 2, 1 }, { 0, 0, 0 }, { -1, -2, -1 } };
if (!img)return;
kernel = create_gaussian_kernel(sigma);
if (!kernel)return;
gaussian_blur(img, kernel, 1);
sgx = sobel_grad(img, sobelx);
gaussian_blur(sgx, kernel, 0);
sgy = sobel_grad(img, sobely);
gaussian_blur(sgy, kernel, 0);
calculate_angle_and_mag(sgx, sgy, &dir, &mag);
non_maximum_supression(dir, mag, &nms, &mag_max);
mag_max = get_real_2d_max(nms);
normalize_real_2d(nms, mag_max);
mag_max = get_real_2d_max(nms);
high_thred = mag_max * high_thred_ratio;
low_thred = high_thred * low_thred_ratio;
int_vec* srv = NULL;
int_vec* scv = NULL;
int_vec* wrv = NULL;
int_vec* wcv = NULL;
int scount = 0;
int wcount = 0;
double_thresholding(nms, low_thred, high_thred, &srv, &scv, &scount, &wrv, &wcv, &wcount);
//edge_track(nms, srv, scv, scount);
edge_track2(nms, srv, scv, scount/*, wcount*/);
cleanup_weak_edge(nms, wrv, wcv, wcount);
normalized_to_unnormalize(nms);
e = real_to_uint8(nms);
destroy_real_2d(&sgx);
destroy_real_2d(&sgy);
destroy_real_2d(&mag);
destroy_real_2d(&nms);
destroy_int_vec(&srv);
destroy_int_vec(&scv);
destroy_int_vec(&wrv);
destroy_int_vec(&wcv);
if (!e)return;
*edge = e->data;
e->data = NULL;
destroy_uint8_2d(&e);
}
三、优化上的考虑:
1、双阈值选取;
2、幅值、梯度、梯度方向可以在一个方法内完成;
3、非最大值抑制可以采用插值实现。
4、高斯模糊算子可以采用先梯度运算,这样可以对一些弱边缘效果很好,matlab就是采用这种实现,但是matlab对后面的实现封装了,没有相应的源码可以学习。
参考资料:
1、https://en.wikipedia.org/wiki/Gaussian_function
2、https://en.wikipedia.org/wiki/Sobel_operator
3、Digital Image Processing 4th Edition Rafael C. Gonzalez.pdf,911-916
特别感谢http://justin-liang.com/tutorials/canny/链接中提供的方法,找过好几个matlab实现代码,效果均没有这里给出的好。