matlab实现灰度图像边缘提取

实验一：编程实现灰度图像空间域边缘提取，至少包括： 梯度算子、Roberts算子、Sobel算子、拉普拉斯算子及LOG算子，并比较不同算法提取边缘的效果及影响因素

代码：

f=rgb2gray(imread('D:\图片\image\liuhangyu.jpg')); %更换你的图片地址
T=20;%阈值
[m,n]=size(f);
%------梯度法-------
f_g=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        f_g(i,j)=abs(f(i+1,j)-f(i,j))+abs(f(i,j+1)-f(i,j));
        if f_g(i,j)<T
            f_g(i,j)=0;
        else
            f_g(i,j)=255;
        end
    end
end
figure(1);
subplot(2,4,1);imshow(uint8(f_g));title('梯度法');
 
%------roberts算子-------
f_r=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        f_r(i,j)=abs(f(i+1,j+1)-f(i,j))+abs(f(i,j+1)-f(i+1,j));
        if f_r(i,j)<T
            f_r(i,j)=0;
        else
            f_r(i,j)=255;
        end
    end
end
%f_r=imbinarize(imfilter(f,r),T);
subplot(2,4,2);imshow(uint8(f_r));title('Roberts算子');
 
%------prewitt算子-------
f_p=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        f_p(i,j)=abs(f(i-1,j-1)+f(i,j-1)+f(i+1,j-1)-f(i-1,j+1)-f(i,j+1)-f(i+1,j+1))+abs(f(i+1,j-1)+f(i+1,j)+f(i+1,j+1)-f(i-1,j-1)-f(i-1,j)-f(i-1,j+1));
        if f_p(i,j)<15
            f_p(i,j)=0;
        else
            f_p(i,j)=255;
        end
    end
end
subplot(2,4,3);imshow(uint8(f_p));title('Prewitt算子');
 
%------sobel算子-------
f_s=zeros(m,n);
for i=2:m-1
    for j=2:n-1
        f_s(i,j)=abs(f(i-1,j-1)+2*f(i,j-1)+f(i+1,j-1)-f(i-1,j+1)-2*f(i,j+1)-f(i+1,j+1))+abs(f(i+1,j-1)+2*f(i+1,j)+f(i+1,j+1)-f(i-1,j-1)-2*f(i-1,j)-f(i-1,j+1));
        if f_s(i,j)<T
            f_s(i,j)=0;
        else
            f_s(i,j)=255;
        end
    end
end
subplot(2,4,4);imshow(uint8(f_s));title('Sobel算子');
 
%------krisch算子-------
k(:,:,1)=[-3,-3,-3;
    -3,0,5;
    -3,5,5];
k(:,:,2)=[-3,-3,5;
    -3,0,5;
    -3,-3,5];
k(:,:,3)=[-3,5,5;
    -3,0,5;
    -3,-3,-3];
k(:,:,4)=[-3,-3,-3;
    -3,0,-3;
    5,5,5];
k(:,:,5)=[5,5,5;
    -3,0,-3;
    -3,-3,-3];
k(:,:,6)=[-3,-3,-3;
    5,0,-3;
    5,5,-3];
k(:,:,7)=[5,-3,-3;
    5,0,-3;
    5,-3,-3];
k(:,:,8)=[5,5,-3;
    5,0,-3;
    -3,-3,-3];
kk=zeros(size(f));
I=double(f);
for i=1:8
    f_k(:,:,i)=conv2(I,k(:,:,i),'same');
    kk=max(kk,f_k(:,:,i));
end
f_kk=imbinarize(kk,600);
subplot(2,4,5);imshow(f_kk);title('Krisch算子');
 
%------LoG算子-------
log1=[0 0 -1 0 0;
    0 -1 -2 -1 0;
    -1 -2 16 -2 -1;
    0 -1 -2 -1 0;
    0 0 -1 0 0];
f_l=conv2(f,log1,'same');
f_ll=imbinarize(abs(f_l),300);
subplot(2,4,6);imshow(f_ll);title('LoG算子');

%------拉普拉斯算子------- 
I=f;
I=im2double(I);
[M,N]=size(I);
B=zeros(size(I));
for x=2:M-1
    for y=2:N-1
        B(x,y)=I(x+1,y+1)+I(x-1,y-1)+I(x-1,y+1)+I(x+1,y-1)+I(x+1,y)+I(x-1,y)+I(x,y+1)+I(x,y-1)-8*I(x,y);
    end
end
I=im2uint8(I);
B=im2uint8(B);
subplot(2,4,7);imshow(I);title('原图');
subplot(2,4,8);imshow(B);title('拉普拉斯算子后的图');

现象：

分析：

通过比较提取边缘效果，可以发现Roberts算子的边缘定位较准确，但是抗噪声能力较弱。Sobel算子有平滑差分作用，它利用像素临近区域的梯度值来计算一个像素的梯度，然后进行取舍。拉普拉斯高斯算子通过对图像进行微分操作实现边缘检测，所以对离散点和噪声较敏感，可以强化突变，并且保留部分物体内部。

实验二：编程实现灰度图像阈值分割，至少包括： 固定阈值分割、最大类间方差（Otsu）方法、最佳熵自动门限方法，并比较不同算法分割的效果及影响因素。

代码

由于算法复杂，此程序需要执行很长时间如果遇到没反应请等待一会即可

tip:担心报错请用下面图片处理

I=imread('D:\图片\image\liuhangyu.jpg');
imshow(I);
%输出直方图
figure;imhist(I);
%人工选定阈值进行分割，选择阈值为120
[width,height]=size(I);
T1=120;
for i=1:width
    for j=1:height
        if(I(i,j)<T1)
            BW1(i,j)=0;
        else 
            BW1(i,j)=1;
        end
    end
end
figure;imshow(BW1),title('人工阈值进行分割');
%自动选择阈值
T2=graythresh(I);
BW2=im2bw(I,T2);%Otus阈值进行分割
figure;imshow(BW2),title('Otus阈值进行分割');

count = imhist(I);                %图像的直方图
[m,n] = size(I);
N = m*n;
L = 256;
count = count/N;                  %每一个像素的分布概率
 
for i = 1:L
    if count(i) ~= 0
        st = i-1;
        break;
    end
end
 
for i = L:-1:1
    if count(i) ~= 0
        nd = i-1;
        break;
    end
end
 
f = count(st+1:nd+1);             %f是每个灰度出现的概率
size(f);
E=[];
for Th = st:nd+1
    Hbt = 0;
    Hwt = 0;
    Pth = sum(count(1:Th+1));
    
    for i = 0:Th            %计算图像背景的熵
        Hbt = Hbt-count(i+1)/Pth*log(count(i+1)/Pth+0.01);  
    end
   
    for i = Th+1:L-1;       %计算图像目标的熵
        Hwt = Hwt-count(i+1)/(1-Pth)*log(count(i+1)/(1-Pth)+0.01);
    end
    E(Th-st+1) = Hbt+Hwt;
end
position = find(E==(max(E)));
Ht = st+position-100;
for i = 1:m
    for j = 1:n
        if a(i,j)>Ht
            a(i,j) = 0;            %对图像分割
        else
            a(i,j) = 255;
        end
    end
end
figure,imshow(I),title('最佳熵自动门限方法');

报错，换上面提示的图片处理
实验现象应有5张图，不要出现数组超出