一、概要:
使用Canny边缘检测算法作为例子,介绍图像的平滑方法和边缘检测。
Canny边缘检测算法分为四步:
step1:用高斯滤波器平滑图像;
step2:用一阶偏导的有限差分来计算梯度的幅值和方向;(在横竖两个方向上计算边缘,再求平方和的开方)
step3:对梯度幅值进行非极大值抑制;
step4:用双阈值算法检测和连接边缘。
demo&效果:
原图(lenna.bmp):
高斯滤波后的图像:
初步求边缘后的图像:
非极大值抑制后的图像:
双阈值检测后的图像:
以下对这四步进行详细介绍:
二、图像平滑与高斯滤波
图像平滑的目的是消除或尽量减少噪声的影响,改善图像的质量。在假定加性噪声是随机独立分布的条件下,利用邻域的平均或加权平均可以有效的抑制噪声干扰。
方法是做邻域计算,比如:
高斯滤波:
采用高斯函数作为加权函数。
原因一:二维高斯函数具有旋转对称性,保证滤波时各方向平滑程度相同;
原因二:离中心点越远权值越小。确保边缘细节不被模糊。
设定σ2和n,确定高斯模板权值。如σ2 =2和n=5,整数化和归一化后得:
代码:
1.直接用系统函数的方法
-
img=imread(
'lenna.bmp');
-
imshow(img);
-
[
m n]=size(img);
-
img=
double(img);
-
%高斯滤波
-
w=fspecial(
'gaussian',[
5
5]); %[
5
5 ]是模板尺寸,默认是[
3
3 ] 模板即上文所提的模板
-
img=imfilter(img,w,
'replicate');
-
figure;
-
imshow(uint8(img))
2.自己实现的高斯滤波(使用上面的模板)
-
img=imread(
'lenna.bmp');
-
img=im2double(img);
-
[m n]=size(img);
-
h=zeros(m,n);
-
for i=
1
:m
-
for j=
1
:n
-
if i<
3
|| i>(m-
3)
|| j<
3
|| j>(n-
3)
-
h(i,j)=img(i,j);
-
else %代码长度原因,分行相加。
-
h(i,j)=h(i,j)+img(i-
2,j-
2)+img(i-
2,j+
2)+img(i+
2,j-
2)+img(i+
2,j+
2);
-
h(i,j)=h(i,j)+(img(i-
1,j-
2)+img(i+
1,j-
2)+img(i-
2,j-
1)+img(i+
2,j-
1)+img(i-
2,j+
1)+img(i+
2,j+
1)+img(i-
1,j+
2)+img(i+
1,j+
2))*
2;
-
h(i,j)=h(i,j)+
3*(img(i,j-
2)+img(i,j+
2)+img(i+
2,j)+img(i-
2,j));
-
h(i,j)=h(i,j)+
4*(img(i-
1,j-
1)+img(i-
1,j+
1)+img(i+
1,j-
1)+img(i+
1,j+
1));
-
h(i,j)=h(i,j)+
6*(img(i,j-
1)+img(i,j+
1)+img(i+
1,j)+img(i-
1,j));
-
h(i,j)=h(i,j)+
7*img(i,j);
-
h(i,j)=h(i,j)/
79;
-
end
-
end
-
end
-
imshow(h,[]);
三、初步边缘检测
对高斯平滑后的图像进行sobel边缘检测。这里需要求横的和竖的还有联合的,所以一共三个需要sobel边缘检测图像。
代码:
-
img=imread(
'lenna.bmp');
-
[m n]=size(img);
-
img=
double(img);
-
w=fspecial(
'gaussian',[5 5]);
-
img=imfilter(img,w,
'replicate');
-
figure;
-
%%sobel边缘检测
-
w=fspecial(
'sobel');
-
img_w=imfilter(img,w,
'replicate'); %求横边缘
-
w=w
'; %转置矩阵
-
img_h=imfilter(img,w,
'replicate'); %求竖边缘
-
img=sqrt(img_w.^
2+img_h.^
2); %平方和再开方 .^表示每一位都和自己乘,不清楚的自己百度
-
figure;
-
imshow(uint8(img))
四、非极大值抑制以及双阈值检测边缘
什么是非极大值抑制?
非极大值抑制是在梯度方向上的极大值。
代码:(觉得看看注释就差不多明白思路了)
-
img=imread(
'lenna.bmp');
-
imshow(img);
-
[m n]=size(img);
-
img=
double(img);
-
%%canny边缘检测的前两步相对不复杂,所以我就直接调用系统函数了
-
%%高斯滤波
-
w=fspecial(
'gaussian',[5 5]);
-
img=imfilter(img,w,
'replicate');
-
figure;
-
imshow(uint8(img))
-
%%sobel边缘检测
-
w=fspecial(
'sobel');
-
img_w=imfilter(img,w,
'replicate'); %求横边缘
-
w=w
';
-
img_h=imfilter(img,w,
'replicate'); %求竖边缘
-
img=sqrt(img_w.^
2+img_h.^
2); %注意这里不是简单的求平均,而是平方和在开方。我曾经好长一段时间都搞错了
-
figure;
-
imshow(uint8(img))
-
%%下面是非极大抑制
-
new_edge=zeros(m,n);
-
for i=
2:m
-1
-
for j=
2:n
-1
-
Mx=img_w(i,j);
-
My=img_h(i,j);
-
-
if My~=
0
-
o=atan(Mx/My); %边缘的法线弧度
-
elseif My==
0 && Mx>
0
-
o=pi/
2;
-
else
-
o=-pi/
2;
-
end
-
-
%Mx处用My和img进行插值
-
adds=get_coords(o);%边缘像素法线一侧求得的两点坐标,插值需要
-
M1=My*img(i+adds(
2),j+adds(
1))+(Mx-My)*img(i+adds(
4),j+adds(
3)); %插值后得到的像素,用此像素和当前像素比较
-
adds=get_coords(o+pi);%边缘法线另一侧求得的两点坐标,插值需要
-
M2=My*img(i+adds(
2),j+adds(
1))+(Mx-My)*img(i+adds(
4),j+adds(
3));%另一侧插值得到的像素,同样和当前像素比较
-
-
isbigger=(Mx*img(i,j)>M1)*(Mx*img(i,j)>=M2)+(Mx*img(i,j)<M1)*(Mx*img(i,j)<=M2); %如果当前点比两边点都大置
1
-
-
if isbigger
-
new_edge(i,j)=img(i,j);
-
end
-
end
-
end
-
figure;
-
imshow(uint8(new_edge))
-
%%下面是滞后阈值处理
-
up=
120; %上阈值
-
low=
100; %下阈值
-
set(
0,
'RecursionLimit',10000); %设置最大递归深度
-
for i=
1:m
-
for j=
1:n
-
if new_edge(i,j)>up &&new_edge(i,j)~=
255 %判断上阈值
-
new_edge(i,j)=
255;
-
new_edge=connect(new_edge,i,j,low);
-
end
-
end
-
end
-
figure;
-
imshow(new_edge==
255)
-
function re=get_coords(angle) %angle是边缘法线角度,返回法线前后两点
-
sigma=
0.000000001;
-
x1=
ceil(
cos(angle+pi/
8)*
sqrt(
2)
-0.5-sigma);
-
y1=
ceil(-
sin(angle-pi/
8)*
sqrt(
2)
-0.5-sigma);
-
x2=
ceil(
cos(angle-pi/
8)*
sqrt(
2)
-0.5-sigma);
-
y2=
ceil(-
sin(angle-pi/
8)*
sqrt(
2)
-0.5-sigma);
-
re=[x1 y1 x2 y2];
-
-
end
-
function nedge=connect(nedge,y,x,low) %种子定位后的连通分析
-
neighbour
=[-1 -1;-
1
0;-
1
1;
0 -
1;
0
1;
1 -
1;
1
0;
1
1]; %八连通搜寻
-
[m n]=size(nedge);
-
for k=
1:
8
-
yy=y+neighbour(k,
1);
-
xx=x+neighbour(k,
2);
-
if yy>=
1 &&yy<=m &&xx>=
1 && xx<=n
-
if nedge(yy,xx)>=low && nedge(yy,xx)~=
255 %判断下阈值
-
nedge(yy,xx)=
255;
-
nedge=connect(nedge,yy,xx,low);
-
end
-
end
-
end
-
-
end
参考资料:
1.http://www.cnblogs.com/tiandsp/archive/2012/12/13/2817240.html
2.南京大学软件学院数字图像处理课程PPT
一些更详细的原理性的东西可以参考:
http://blog.csdn.net/likezhaobin/article/details/6892176