（亲测可用）基于matlab的用自写函数来实现图像的灰度处理sobel canny算子边缘检测

I=imread(‘E:\matlab\toolbox\images\imdata\tupian4.jpg’); %图片的位置是自己设置的

subplot(2,5,1);
imshow(I);
title(‘原始图像’);

G=I(:,:,2);
subplot(2,5,2);
imshow(G);
title(‘灰度处理后图像’); %彩色图像需要灰度处理，黑白图像则不需要

c=avefilt(G,3); %对G用3x3的模板进行滤波
subplot(2,5,3);
imshow(uint8( c ));
title(‘均值滤波后的图像’);

hx=[-1 -2 -1;0 0 0 ;1 2 1]; %生产sobel垂直梯度模板，这里用的是自定义的模板
hy=hx’; %生产sobel水平梯度模板

gradx=filter2(hx,c,‘same’);
gradx=abs(gradx); %计算图像的sobel垂直梯度

grady=filter2(hy,c,‘same’);
grady=abs(grady); %计算图像的sobel水平梯度

grad=gradx+grady;                  %得到图像的sobel梯度
subplot(2,5,4);
imshow(grad,[ ]);
title('图像的sobel梯度');

%canny算子--------------------------------------------------------
[m,n]=size(G);
src_gaussion=my_gaussian(G,3,9); %高斯滤波

w=fspecial(‘sobel’); %这里的sobel算子是与上面不同的算法
src_sobel_x=imfilter(src_gaussion,w,‘replicate’); %求横边缘
w=w’;
src_sobel_y=imfilter(src_gaussion,w,‘replicate’); %求竖边缘
src_sobel_x=double(src_sobel_x);
src_sobel_y=double(src_sobel_y);
src_sobel=sqrt(src_sobel_x.^{2+src_sobel_y.}2); %注意平方和在开方。
% figure;
subplot(2,5,5);
imshow(uint8(src_sobel_x))

title(‘这是sobel进行梯度处理后的图像’);
new_angel_img=zeros(m,n);
for i=1:m
for j=1:n
Mx=src_sobel_x(i,j);
My=src_sobel_y(i,j);

    if My~=0  
        o=atan(Mx/My);      %边缘的法线弧度  
    elseif My==0 && Mx>0  
        o=pi/2;  
    else  
        o=-pi/2;              
    end  
      
    if( ( o>=( -22.5/180*pi ) && o<( 22.5/180*pi ) ) || ( o>=( 157.5/180*pi ) && o<=pi ) || ( o<=( -157.5/180*pi ) && o>=-pi ) )  
        d1=0;  
    elseif (  ( o>= ( 22.5/180*pi ) && o<( 67.5/180*pi ) ) ||  (   o>=( -157.5/180*pi ) && o<(-112.5/180*pi )  )   )  
        d1=-45;  
    elseif ( ( o>=(67.5/180*pi) && o<(112.5/180*pi) ) || ( o>=(-112.5/180*pi) && o<-67.5/180*pi ) )     
        d1=90;  
    else        
        d1=45;  
     end  
        new_angel_img(i,j)=d1;  
end  

end

canny_dst=src_sobel;
%%%%%%%%%%%%%%%%%%非最大值抑制%%%%%%%%%%%%%%%%%%%%%%%5
for i=3:m-2
for j=3:n-2
if( new_angel_img(i,j)0 )
A=[canny_dst(i-1,j),canny_dst(i+1,j)];
if( abs(canny_dst(i,j))< max( abs(A) ) )
canny_dst(i,j)=0;
end
elseif( new_angel_img(i,j)-45 )
A=[canny_dst(i-1,j-1),canny_dst(i+1,j+1)];
if( abs(canny_dst(i,j))< max( abs(A) ))
canny_dst(i,j)=0;
end
elseif( new_angel_img(i,j)==90 )
A=[canny_dst(i,j-1),canny_dst(i,j+1)];
if( abs(canny_dst(i,j))< max( abs(A) ) )
canny_dst(i,j)=0;
end
elseif( new_angel_img(i,j)==45 )
A=[canny_dst(i-1,j+1),canny_dst(i+1,j-1)];
if( abs(canny_dst(i,j))< max( abs(A) ) )
canny_dst(i,j)=0;
end
end
end
end

% figure(2);
subplot(2,5,6);
imshow(canny_dst);

title(‘这是非最大值抑制11’);
gh=zeros(m,n);
gl=zeros(m,n);
gl1=zeros(m,n);
for i=1:m
for j=1:n
if( canny_dst(i,j)>=100 )
gh(i,j)=1;
else
gh(i,j)=0;
end
if( canny_dst(i,j)>=30)
gl(i,j)=1;
else
gl(i,j)=0;
end
end
end
gl1=gl;
gl1=logical(gl1);
gh=logical(gh);
% figure(3);
subplot(2,5,7);
imshow(gl1);
title(‘这是非最大值抑制22’);

% figure(4);
subplot(2,5,8);
imshow(gh);
title(‘这是递归连接后的33’);
%%%%%%%%%%%%%%%%%%%%%%%%递归连接%%%%%%%%%%%%%%%%%%%
up=100; %上阈值
low=33; %下阈值
set(0,‘RecursionLimit’,10000); %设置最大递归深度
for i=1:m
for j=1:n
if canny_dst(i,j)>up &&canny_dst(i,j)~=255 %判断上阈值
canny_dst(i,j)=255;
canny_dst=connect(canny_dst,i,j,low);
end
end
end
% figure(5);
subplot(2,5,9);
imshow(canny_dst==255);

title(‘这是递归连接后的’);
BW2 = edge(src_gaussion,‘canny’);
% figure(6);
subplot(2,5,10);
imshow(BW2);
title(‘这是递归连接后的2’);

%所需函数————————————————————————————————————

function d=my_gaussian(src,n,k)

n1=(n+1)/2; %%%%%%%%%%%%%计算高斯滤波模版中心
[m,l]=size(src);
b=zeros(n,n); %%%%%%%%%%%%%空模版
I=zeros(m,l);
d=zeros(m,l);
img=zeros(m,l);
%%%%%%%%%%%%%%%%%计算高斯滤波模版的各个系数
for i=1:n
for j=1:n
b(i,j)=(exp(-((i-n1)^2 +(j-n1)^2)/(2k)))/(2pi*k);
end
end

b = b/sum(b( : )); %%%%%%%%%%%%%%%%%%%%%%对系数进行归一化 （这个地方必须归一化，不然整体灰度会底） 删除括号里面的空格
I=double(src);
img=conv2(I,b,‘same’); %%%%%%%%%%%%图像与高斯滤波模版卷积，输出与输入图像大小相同的图像
d=uint8(img); %%%%%%%%%%转为无符号整型，用来显示

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

function d=avefilt(x,n)
a(1:n,1:n)=1; %a即n×n模板,元素全是1
p=size(x); %输入图像是p×q的,且p>n,q>n
x1=double(x);
x2=x1;
%A(a:b,c:d)表示A矩阵的第a到b行,第c到d列的所有元素
for i=1:p(1)-n+1
for j=1:p(2)-n+1
c=x1(i:i+(n-1),j:j+(n-1)).a; %取出x1中从(i,j)开始的n行n列元素与模板相乘
s=sum(sum©); %求c矩阵(即模板)中各元素之和
x2(i+(n-1)/2,j+(n-1)/2)=s/(nn); %将模板各元素的均值赋给模板中心位置的元素
end
end
% %未被赋值的元素取原值
d=uint8(x2);
end