基于梯度的边缘检测算子(sobel、canny、roberts、prewitt)原理及MATLAB代码

最新推荐文章于 2024-04-23 23:40:07 发布

weixin_45858794

最新推荐文章于 2024-04-23 23:40:07 发布

阅读量7.2k

点赞数 7

分类专栏：图像处理文章标签：边缘检测计算机视觉

本文链接：https://blog.csdn.net/weixin_45858794/article/details/106159560

版权

图像处理专栏收录该内容

2 篇文章 0 订阅

订阅专栏

一、Sobel算子
sobel算子的核心是像素矩阵的卷积，而卷积本质就是对指定的图像区域的像素值进行加权求和的过程，其计算过程为图像区域中的每个像素值分别与卷积模板的每个元素对应相乘，将卷积的结果作求和运算，运算到的和就是卷积运算的结果。
sobel算子包含垂直（左）和水平（右）两个方向的卷积模板在这里插入图片描述水平梯度模板在这里插入图片描述
若A为原始图像，则

改变后的灰度值有两种计算方式：

最后设置一个阈值，运算后的像素值大于该阈值输出为0，小于该阈值输出为255。

sobel算子matlab代码：

 %RGB_YCbCr 
  clc; 
  clear all; 
  close all; 
  RGB_data = imread('微信图片_20191227124400.jpg');% 
  R_data =    RGB_data(:,:,1); 
  G_data =    RGB_data(:,:,2);
  B_data =    RGB_data(:,:,3);
  imshow(RGB_data);
  title('原图像');
 [ROW,COL, DIM] = size(RGB_data); 
 Y_data = zeros(ROW,COL);
 Cb_data = zeros(ROW,COL);
 Cr_data = zeros(ROW,COL);
 Gray_data = RGB_data; 
 for r = 1:ROW 
    for c = 1:COL
         Y_data(r, c) = 0.299*R_data(r, c) + 0.587*G_data(r, c) + 0.114*B_data(r, c);%彩色图像变为灰度图像公式Gray = R*0.299 + G*0.587 + B*0.114
         Cb_data(r, c) = -0.172*R_data(r, c) - 0.339*G_data(r, c) + 0.511*B_data(r, c) + 128;
        Cr_data(r, c) = 0.511*R_data(r, c) - 0.428*G_data(r, c) - 0.083*B_data(r, c) + 128;
     end
 end 
 Gray_data(:,:,1)=Y_data;
 Gray_data(:,:,2)=Y_data;
 Gray_data(:,:,3)=Y_data; 
 figure;
 imshow(Gray_data);
 title('灰度图像');
 %Median Filter
 imgn = imnoise(Gray_data,'salt & pepper',0.02); %当噪声类型是’salt & pepper’的时候，第三个参数的意思是噪声密度，比如0.1，那么总像素个数的10%为黑白点，当然是黑点还是白点都是随机的。
 figure;
 imshow(imgn);
  title('椒盐噪声图像');
 Median_Img = Gray_data;
 for r = 2:ROW-1
     for c = 2:COL-1
         median3x3 =[imgn(r-1,c-1)    imgn(r-1,c) imgn(r-1,c+1)
                     imgn(r,c-1)      imgn(r,c)      imgn(r,c+1)
                     imgn(r+1,c-1)      imgn(r+1,c) imgn(r+1,c+1)];
         sort1 = sort(median3x3, 2, 'descend');%对矩阵的每行进行降序排列
         sort2 = sort([sort1(1), sort1(4), sort1(7)], 'descend');%对矩阵第一列进行降序排列
         sort3 = sort([sort1(2), sort1(5), sort1(8)], 'descend');%对矩阵第二列进行降序排列
         sort4 = sort([sort1(3), sort1(6), sort1(9)], 'descend');%对矩阵第三列进行降序排列
         mid_num = sort([sort2(3), sort3(2), sort4(1)], 'descend');%对第一列最小值，第二列中间值，第三列最大值进行降序排列
         Median_Img(r,c) = mid_num(2);%取中间值
     end
 end
 figure;
 imshow(Median_Img);
  title('中值滤波后图像');
 %Sobel_Edge_Detect
 Median_Img = double(Median_Img);
 Sobel_Threshold = 160;
 Sobel_Img = zeros(ROW,COL);
 for r = 2:ROW-1
     for c = 2:COL-1
         Sobel_x = Median_Img(r-1,c+1) + 2*Median_Img(r,c+1) + Median_Img(r+1,c+1) - Median_Img(r-1,c-1) - 2*Median_Img(r,c-1) - Median_Img(r+1,c-1);
         Sobel_y = Median_Img(r-1,c-1) + 2*Median_Img(r-1,c) + Median_Img(r-1,c+1) - Median_Img(r+1,c-1) - 2*Median_Img(r+1,c) - Median_Img(r+1,c+1);
         Sobel_Num = abs(Sobel_x) + abs(Sobel_y);
         %Sobel_Num = sqrt(Sobel_x^2 + Sobel_y^2);
         if(Sobel_Num > Sobel_Threshold)            
             Sobel_Img(r,c)=0;
         else
             Sobel_Img(r,c)=255;
         end
     end
 end 
 figure;
 imshow(Sobel_Img);
  title('Sobel边缘检测图像');

Sobel算法只采用水平和垂直2个方向模板，对其他方向的边缘不敏感，对于纹理较复杂、斜向边缘较多的图像轮廓提取不是很理想，抗噪声能力也比较低。

改进方法：
在原有水平和垂直方向上增加了两个对角线方向上的模板

二、Canny 算子
Canny 算子进行边缘检测步骤：
1）对图像进行二维高斯滤波；
通过使用高斯滤波器与图像进行卷积，来减少噪声对边缘检测结果的影响和平滑图像，以减少边缘检测器上明显的噪声影响。
2）通过一阶微分计算图像的灰度梯度幅值和方向；
通过梯度方向确定边缘的方向，就可以把边缘的梯度方向大略分成几种角度(如0°、45°、90°和135 ), 并可以找到这个像素梯度方向的邻接像素。该梯度值的计算和sobel计算过程一样，改变后的灰度值有两种计算方式：
在这里插入图片描述
3）对计算出的梯度幅值进行非极大值抑制（Non-MaximaSuppression, NMS）；
由于梯度值最大不一定为边缘点，需要对一阶微分计算后的图像数据进行非极大值抑制，只保留幅值局部变化最大的点。将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较。如果当前像素的梯度强度与另外两个像素相比最大，则该像素点保留为边缘点，否则该像素点将被抑制。
4）通过人为设定的高低阈值确定图像的边缘。
为了解决噪声和颜色变化引起的一些边缘像素，必须用弱梯度值过滤边缘像素，并保留具有高梯度值的边缘像素，可以通过选择高低阈值来实现。如果边缘像素的梯度值高于高阈值，则将其标记为强边缘像素；如果边缘像素的梯度值小于高阈值并且大于低阈值，则将其标记为弱边缘像素；如果边缘像素的梯度值小于低阈值，则会被抑制。阈值的选择取决于给定输入图像的内容。

canny算子matlab代码：

clc;
clear all;
image = imread('微信图片_20191227124400.jpg'); 
image = rgb2gray(image);
subplot(221); 
imshow(image); 
title('原始图像'); 

image = double(image)/256; 
[m,n] = size(image); 
w = fspecial('gaussian');
image_1 = imfilter(image,w,'replicate');%滤波
subplot(222);
imshow(int8(256*image_1)); 
title('高斯滤波后的图像');  % 梯度计算

op = fspecial('sobel')/4;  % 用Sobel算子来求导数 
x = op'; y =op; 
bx = imfilter(image_1,x,'replicate');
by = imfilter(image_1,y,'replicate'); 
b_abs = sqrt(bx.*bx+by.*by);        % 求梯度的幅值 
b_angle = angle(by-i*bx);
b_ang = b_angle/3.1416*180;         % 求梯度的方向 

% 梯度方向确定
for r = 1:m   
    for c = 1:n        
        if((b_ang(r,c)>=22.5 && b_ang(r,c)<67.5) || (b_ang(r,c)>=-157.5 && b_ang(r,c)<-112.5))           
                dir(r,c) = 1;        
        else if ((b_ang(r,c)>=67.5 && b_ang(r,c)<112.5)|| (b_ang(r,c)>=-112.5 && b_ang(r,c)<-67.5))        
                dir(r,c) = 2;        
            else if ((b_ang(r,c)>=112.5 && b_ang(r,c)<157.5)|| (b_ang(r,c)>=-67.5 && b_ang(r,c)<-22.5))           
                    dir(r,c) = 3;        
                else
                    dir(r,c) = 0;  
                end
                end
        end
    end
end

% 遍历图像
        b_ab = [zeros(m,1),b_abs,zeros(m,1)];    % 串联矩阵
        b_ab = [zeros(1,n+2);b_ab;zeros(1,n+2)]; %相当于在求梯度后的像素值加了一周0
        for r = 2:m+1     
            for c = 2:n+1        
                switch dir(r-1,c-1)             
                    case 0                 
                        if((b_ab(r,c)<b_ab(r+1,c))| (b_ab(r,c)<b_ab(r-1,c)))                   
                            b1(r-1,c-1) = 0;                
                        else
                            b1(r-1,c-1) = b_ab(r,c);                
                        end
                    case 1               
                        if((b_ab(r,c)<b_ab(r+1,c-1))| (b_ab(r,c)<b_ab(r-1,c+1)))                
                            b1(r-1,c-1) = 0;                
                        else
                            b1(r-1,c-1) = b_ab(r,c);             
                        end
                    case 2               
                        if((b_ab(r,c)<b_ab(r,c-1))| (b_ab(r,c)<b_ab(r,c+1)))                 
                            b1(r-1,c-1) = 0;              
                        else
                            b1(r-1,c-1) = b_ab(r,c);             
                        end
                    case 3             
                        if((b_ab(r,c)<b_ab(r-1,c-1))| (b_ab(r,c)<b_ab(r+1,c+1)))              
                            b1(r-1,c-1) = 0;               
                        else
                            b1(r-1,c-1) = b_ab(r,c);                
                        end
                end
            end
        end
        
        for r = 1:m     
            for c = 1:n        
                if (b1(r,c)>0.24)           
                    b2(r,c) = 1;       
                else
                    b2(r,c) = 0;       
                end
            end
        end
        
        for r = 1:m    
            for c = 1:n       
                if(b1(r,c)>0.36)           
                    b3(r,c)=1;       
                else
                    b3(r,c) = 0;       
                end
            end
        end
        
        image_2 = b3; 
        for k = 1:10    
            for r = 2:m-1     
                for c = 2: n-1          
                    if (b2(r,c)==1 && (image_2(r,c)==1||image_2(r+1,c)==1 ...
                        ||image_2(r+1,c-1)==1||image_2(r+1,c+1)==1||image_2(r,c-1)==1 ...
                    ||image_2(r,c+1)==1||image_2(r-1,c-1)==1||image_2(r-1,c)==1 ...
                    ||image_2(r-1,c+1)==1))                 
                        image_2(r,c) = 1;        
                    else
                        image_2(r,c) = 0;          
                    end
                end
            end
        end
        subplot(223);
        imshow(image_2);
        title('Canny算子检测后的图像');

由于梯度计算对噪声非常敏感，并且采用高斯滤波需要人为设定方差，所以容易出现漏检现象。

改进方法：
利用双线性滤波代替高斯滤波，同时在非极大值抑制时将梯度幅值极大值点进行保留并作为候选边缘点，并且进行自适应阈值选取。


clear all;
close all;
%%%%%%%%%%%%%%载入图像%%%%%%%%%%%%%%
X = zeros(200);                                       % 定义一个灰度块矩阵A
X(1:100,1:100) = 85;
X(101:200,101:200) = 170;
X(101:200,1:100) = 255;
[dy, dx]=gradient(X,.1,.1);                           % 求梯度
figure('NumberTitle', 'off', 'Name', '原灰度块1');
imshow(uint8(X))  % 显示原灰度块1
[Ix,Iy]=size(X);           %图像长宽
%%%%%%%%%%%%%%%%%%双边滤波器%%%%%%%%%%%%%%%%%%%%%%%
r=15;        %模板半径  
imgn=zeros(Ix+2*r+1,Iy+2*r+1);  
imgn(r+1:Ix+r,r+1:Iy+r)=X;  
imgn(1:r,r+1:Iy+r)=X(1:r,1:Iy);                 %扩展上边界  
imgn(1:Ix+r,Iy+r+1:Iy+2*r+1)=imgn(1:Ix+r,Iy:Iy+r);    %扩展右边界
imgn(Ix+r+1:Ix+2*r+1,r+1:Iy+2*r+1)=imgn(Ix:Ix+r,r+1:Iy+2*r+1);    %扩展下边界
imgn(1:Ix+2*r+1,1:r)=imgn(1:Ix+2*r+1,r+1:2*r);       %扩展左边界  
sigma_d=7; sigma_r=0.2;  
[x,y] = meshgrid(-r:r,-r:r); 
w1=exp(-(x.^2+y.^2)/(2*sigma_d^2));     %以距离作为自变量高斯滤波器  
for i=r+1:Ix+r    
    for j=r+1:Iy+r                
        w2=exp(-(imgn(i-r:i+r,j-r:j+r)-imgn(i,j)).^2/(2*sigma_r^2)); %以周围和当前像素灰度差值作为自变量的高斯滤波器       
        w=w1.*w2;        
        s=imgn(i-r:i+r,j-r:j+r).*w;     
        imgn(i,j)=sum(sum(s))/sum(sum(w));    
    end
end
figure;
imshow(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)));
title('双边滤波器图像');%显示滤波后的图像

%%%%%%%%%%%%%%二级尺度%%%%%%%%%%%%%%    
m=1.0;
delta=2^m;
%%%%%%%%%高斯函数X,Y方向的偏导%%%%%%%
N=20;                     %长度和幅度
A=-1/sqrt(2*pi);
phi_x=zeros(N,N);
phi_y=zeros(N,N);
for i=1:N                 %偏导
    for j=1:N
        x=i-(N+1)/2;
        y=j-(N+1)/2;
        phi_x(i,j)=A*(x/delta^2).*exp(-(x*x+y*y)/(2*delta^2));
        phi_y(i,j)=A*(y/delta^2).*exp(-(x*x+y*y)/(2*delta^2));
    end
end
%%%%%%%%%图象行、列卷积%%%%%%%%%%%%%%
 W_x=conv2(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)),phi_x,'same');
W_y=conv2(mat2gray(imgn(r+1:Ix+r,r+1:Iy+r)),phi_y,'same');
%%%%%%%%%%%%求梯度%%%%%%%%%%%%%%%%%%%
Mod_G=zeros(Ix,Iy);
for i2=1:Ix
    for j2=1:Iy
        Mod_G(i2,j2)=sqrt(W_x(i2,j2)*W_x(i2,j2)+W_y(i2,j2)*W_y(i2,j2));
    end
end
M=fix(Mod_G);
%%%%%%%%%%%%求幅角%%%%%%%%%%%%%%%%%%%
Arg=zeros(Ix,Iy);
for i3=1:Ix
    for j3=1:Iy
        angle=atan(W_y(i3,j3)/W_x(i3,j3))*180/pi;      %求反正切
        if ((W_x(i3,j3)>0)&&(W_y(i3,j3)>0))            %第一象限
            angle=angle+0;
        elseif ((W_x(i3,j3)<0)&&(W_y(i3,j3)>0))        %第二象限
            angle=angle+180;
        elseif ((W_x(i3,j3)<0)&&(W_y(i3,j3)<0))        %第三象限
            angle=angle+180;
        elseif ((W_x(i3,j3)>0)&&(W_y(i3,j3)<0))        %第四象限
            angle=angle+360;                       
        elseif (abs(W_y(i3,j3)/W_x(i3,j3))>100000000)  %大值90度
            angle=90;
        elseif (abs(W_y(i3,j3)/W_x(i3,j3))==0)         %零值0度
            angle=0;
        end
        if(angle>180)                                  %幅角限制在0-180                        
            angle=angle-180;
        end
        Arg(i3,j3)=angle;
    end
end
%%%%%%%%%%%%图像边界%%%%%%%%%%%%%%%%%%%%
%%%%%%沿着梯度方向检测小波变换系数模的局部极大值点即可得到图像的边缘点%%%%%%%%%
Edge=zeros(Ix,Iy);
for i4=2:(Ix-1)
    for j4=2:(Iy-1)
        if ((Arg(i4,j4)>=0&&Arg(i4,j4)<=22.5))
            if (Mod_G(i4,j4)>Mod_G(i4+1,j4)&&Mod_G(i4,j4)>Mod_G(i4-1,j4))%水平
                Edge(i4,j4)=Mod_G(i4,j4);
            end
        elseif ((Arg(i4,j4)>=(90-22.5)&&Arg(i4,j4)<=(90+22.5)))
            if (Mod_G(i4,j4)>Mod_G(i4,j4+1)&&Mod_G(i4,j4)>Mod_G(i4,j4-1))%45度
                Edge(i4,j4)=Mod_G(i4,j4);
            end
        elseif ((Arg(i4,j4)>=(45-22.5)&&Arg(i4,j4)<=(45+22.5)))
            if (Mod_G(i4,j4)>Mod_G(i4+1,j4+1)&&Mod_G(i4,j4)>Mod_G(i4-1,j4-1))%垂直
                Edge(i4,j4)=Mod_G(i4,j4);
            end
        else  
            if (Mod_G(i4,j4)>Mod_G(i4+1,j4-1)&&Mod_G(i4,j4)>Mod_G(i4-1,j4+1))%135度
                Edge(i4,j4)=Mod_G(i4,j4);
            end
        end
    end
end
MAX_E=max(Edge(:));
if (MAX_E>0)
    Edge=Edge/MAX_E;                           %归一化
end
%%%%%%%%%%%根据直方图获取高低阈值%%%%%%%%%%
count=0;                                       %累计直方图统计可能是边缘点的个数
hist=zeros(1,1024);
for k1=1:Ix
    for k2=1:Iy
        if (Edge(k1,k2)~=0)
            hist(1,M(k1,k2)+1)=hist(1,M(k1,k2)+1)+1;
            count=count+1;
        end
    end
end
for k3=1:1024
    if hist(1,k3)~=0
        nmaxmag=k3;
    end
end
p=0.7;
dot=ceil(p*count);
k4=1;
nedgenb=hist(1,1);
while (k4<(nmaxmag-1)&&(nedgenb<dot))
    k4=k4+1;
    nedgenb=nedgenb+hist(1,k4);
end
k4=k4/1000;
high=k4;                                       %高阈值
low=0.4*high;                                  %低阈值,一般高为低的2.5倍
%%%%%%%%%%%%边缘点判定%%%%%%%%%%%%%%%%%%
for i5=1:Ix
    for j5=1:Iy
        if (Edge(i5,j5)>=high)                 %高于高阈值一定是边缘
            Edge(i5,j5)=255;
        elseif (Edge(i5,j5)<=low)              %低于高阈值一定是边缘
            Edge(i5,j5)=0;
        else                                   %高低阈值之间判断8邻域是否高于阈值否则不是边缘
            if((Edge(i5-1,j5)>high)&&(Edge(i5,j5-1)>high)&&(Edge(i5+1,j5)>=high)&&(Edge(i5,j5+1)>=high)...
                    &&(Edge(i5-1,j5-1)>high)&&(Edge(i5+1,j5-1)>high)&&(Edge(i5-1,j5+1)>high)&&(Edge(i5+1,j5+1)>high))
               Edge(i5,j5)=255;
            else
               Edge(i5,j5)=0;
            end
        end
    end
end

figure;
imshow(Edge);
title('边缘检测图像')

三、Roberts算子
Roberts算子和Sobel算子原理基本一样。只是Roberts算子使用的偶数模板，如2*2，通过对角线上的值对图像像素进行偏移，其垂直（左）和水平（右）卷积模板为：
在这里插入图片描述
若A为原始图像，则

梯度两种计算方式为：

通常为了简化计算，而采用第二种。
最后设置一个阈值，运算后的像素值大于该阈值输出为0，小于该阈值输出为255。

roberts算子matlab代码：

clc;
clear all;
img = imread('微信图片_20191227124400.jpg'); % 读取图像
figure;
imshow(img);
title('原始图像');
img_gray = rgb2gray(img); % 转换成灰度图
figure;
imshow(img_gray);
title('灰度图像');
[m,n] = size(img_gray); % 得到图像的大小
new_img_gray = img_gray; % 新建一个一样大的图像 
pxValue = 0; % roberts计算所得到的像素值 
% 对边界象素操作 
threshold_value=30;     
    for i=1:m-1        
        for j=1:n-1          
            pxValue = abs(img_gray(i,j)-img_gray(i+1,j+1))+...  
                abs(img_gray(i+1,j)-img_gray(i,j+1));      
            if(pxValue > threshold_value)      
                new_img_gray(i,j) = 0;          
            else
                new_img_gray(i,j) =255;      
            end
        end
    end
    figure;
    imshow(new_img_gray);  
    title('roberts边缘检测图像');

roberts算子对图像边缘检测而言定位精度高，在水平和垂直方向效果较好。但该算子对噪声较敏感，无法消除局部干扰，对图像中目标和背景灰度差异表现并不显著的弱边缘却很难检测识别，这将导致提取的目标边缘出现间断。并且该算子的阈值需要人为设定，使得传统的 roberts算子在提取不同目标轮廓时，具有很大的局限性。

改进方法：
采用３×３邻域代替roberts算法中２×２邻域来计算梯度幅值；并利用图像块之间相似性的三维块匹配的去噪模型，提高roberts算子的检测精度和抗噪性能；通过最佳阈值迭代方法代替人为指定阈值来获取最佳分割阈值，有效地提取图中目标轮廓。

四、Prewitt算子
Prewitt算子采用 3*3卷积模板对区域内的像素值进行计算。
其垂直（左）和水平（右）卷积模板为：
在这里插入图片描述
若A为原始图像，则

梯度两种计算方式为：

最后设置一个阈值，运算后的像素值大于该阈值输出为0，小于该阈值输出为255。

prewitt算子matlab代码：

I=imread('微信图片_20191227124400.jpg');
I=im2double(I);
figure;
imshow(I);
title('原始图像'); 
img_gray = rgb2gray(I); % 转换成灰度图
figure;
imshow(img_gray);
title('灰度图像');
[height width R]=size(img_gray); 
threshold_value=0.5;     
for i=2:height-1    
    for j=2:width-1      
        Dx=[I(i+1,j-1)-I(i-1,j-1)]+[I(i+1,j)-I(i-1,j)]+[I(i+1,j+1)-I(i-1,j+1)];      
        Dy=[I(i-1,j+1)-I(i-1,j-1)]+[I(i,j+1)-I(i,j-1)]+[I(i+1,j+1)-I(i+1,j-1)];       
        P(i,j)=sqrt(Dx^2+Dy^2);    
          if (P(i,j)>threshold_value)        
            P(i,j)=0;       
          else
              P(i,j)=255;
    end
    end
end
figure;
imshow(P);
title('prewitt边缘检测图像');

由于Prewitt 算子模板只能对水平和垂直方向的边缘进行检测，无法检测斜侧方向，使边缘检测不完整。

改进方法：
可以通过添加左右斜侧方向的模板，使边缘检测更加完整。

参考文献：[1]https://blog.csdn.net/ahilll/article/details/82050358
[2]https://blog.csdn.net/Jacky_Ponder/article/details/64126411?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.nonecase
[3]https://blog.csdn.net/Prince_IT/article/details/90205259
[4]https://blog.csdn.net/weixin_43788499/article/details/84502478?locationNum=7&fps=1
[5]https://blog.csdn.net/baidu_21578557/article/details/51786249?locationNum=2&fps=1
[6]https://blog.51cto.com/3754839/1401776

weixin_45858794

关注

7
点赞
踩
110

收藏

觉得还不错? 一键收藏
1
评论
基于梯度的边缘检测算子(sobel、canny、roberts、prewitt)原理及MATLAB代码

一、sobel算子sobel算子的核心是像素矩阵的卷积，而卷积本质就是对指定的图像区域的像素值进行加权求和的过程，其计算过程为图像区域中的每个像素值分别与卷积模板的每个元素对应相乘，将卷积的结果作求和运算，运算到的和就是卷积运算的结果。sobel算子包含水平和垂直两个方向的卷积模板...
复制链接

扫一扫