matlab 分别用sobel prewitt roberts laplacian算子进行边缘提取

pic2=imread（‘fig3.tif‘）;
subplot（231）;
imshow（pic2）;
title（‘原始图像‘）;
[ROWCOL]=size（pic2）;                   %取表示图像的矩阵行列值数
count=imhist（pic2）;                     %直方图统计                    %图像矩阵大小
N=ROW*COL;                             %图像像素个数
L=256;                                 %指定图像灰度级为256级
count=count/N;                         %各灰度级出现的概率
%================寻找出现概率不为0的最小灰度值=============
for i=2:L                          
    if count（i）~=0
        st=i-1;
        break;
    end
end
%================寻找出现概率不为0的最大灰度值=============
for i=L:-1:1
    if count（i）~=0
        nd=i-1;
        break;
    end
end
%=========================计算阈值========================
f=count（st+1:nd+1）;
p=st;q=nd-st;                      %p和q分别是灰度的起始和结束值
u=0;
for i=1:q
    u=u+f（i）*（p+i-1）;
    ua（i）=u;
end
for i=1:q
    w1（i）=sum（f（1:i））;
end
d=（u*w1-ua）.^2./（w1.*（1-w1））;       %求出不同i值时类间方差
[ytp]=max（d）;                      %求出最大方差对应的灰度级
th=tp+p                            %threshold
if th>200
    th=100;
end;
% if th<40
%     th=45;
% end;
W=zeros（ROWCOL）;                   %新建一幅ROW×COL大小的图

X=double（pic2）;

%%%从第2行第2列开始到倒数第2行倒数第2列

for i=2:ROW-1

    for j=2:COL-1
        %%%%sobel算子
        %%%%Dx是对X的偏导数，检测水平边缘，梯度方向是垂直的.
        Dx=X（i+1j-1）+2*X（i+1j）+X（i+1j+1）-X（i-1j-1）-2*X（i-1j）-X（i-1j+1）;
        %%%%Dy是对Y的偏导数，检测垂直边缘的梯度方向是水平的.
        Dy=X（i-1j+1）+2*X（ij+1）+X（i+1j+1）-X（i-1j-1）-2*X（ij-1）-X（i+1j-1）;
        %%%梯度的幅值      
        Mag=sqrt（Dx*Dx+Dy*Dy）;
        %%%梯度的方向角Dy有取0的可能.
        %%%a（ij）=round（180*（atan（Dx/Dy）））;
        %%%梯度取整赋值给输出图像Wanted=W;
        W（ij）=round（Mag）;
    end
end

for i=1:ROW
    for j=1:COL
        if W（ij）>th
           W（ij）=255;
        else
           W（ij）=0;
        end
    end
end
subplot（232）;
imshow（W）;
title（‘Sobel算子边缘提取‘）;

for i=2:ROW-1

    for j=2:COL-1
        %%%%prewitt算子
        Dx=X（i+1j-1）+X（i+1j）+X（i+1j+1）-X（i-1j-1）-X（i-1j）-X（i-1j+1）;
        Dx=Dx*4/3;
        Dy=X（i-1j+1）+X（ij+1）+X（i+1j+1）-X（i-1j-1）-X（ij-1）-X（i+1j-1）;  
        Dy=Dy*4/3;
        Mag=sqrt（Dx*Dx+Dy*Dy）;
        %%%梯度的方向角Dy有取0的可能.
        %%%a（ij）=round（180*（atan（Dx/Dy）））;
        %%%梯度取整赋值给输出图像Wanted=W;
        W（ij）=round（Mag）;
    end
end

%%%阈值人为设定随不同图像不同目的改变.
%%%大于阈值认为是边缘设为最亮;小于则为背景最黑

for i=1:ROW
    for j=1:COL
        if W（ij）>th
           W（ij）=255;
        else
           W（ij）=0;
        end
    end
end
subplot（233）;
imshow（W）;
title（‘Prewitt算子边缘提取‘）;

for i=1:ROW-1

    for j=1:COL-1
        %%%%prewitt算子
        Dx=X（ij