图像分割——边缘检测——LoG（Marr-Hildreth算法）（Matlab）_marr-hildreth算子matlab

clc;
clear all;
close all;
%边缘测测试图像（Detection of Edge）
I=im2double(imread('D:\Gray Files\10-16.tif'));
[M,N]=size(I);
 
%%
%=============================边缘检测(四)=================================
% Marr-Hildreth算法,LoG
%-------------------------产生高斯低通滤波器（高斯内核）-------------------
n=25;
n_l=floor(n/2);
Gaussian=zeros(n,n);
delta=4;
for i=-n_l:n_l
    for j=-n_l:n_l
        %计算点（x,y）到中心点的距离
        D=(i)^2+(j)^2;
        %用滤波器乘以主函数 
        Gaussian(i+n_l+1,j+n_l+1)=exp(-D/(2*delta^2));
    end
end
%找到抽样最小值
m=min(Gaussian(:));
Gaussian=floor(Gaussian/m);
%忽略半径大于3*delta的取值，即置零
for i=-n_l:n_l
    for j=-n_l:n_l
        %计算点（x,y）到中心点的距离
        d=sqrt((i)^2+(j)^2);
        %将大于3*delta的取值置零
        if d>3*delta
            Gaussian(i+n_l+1,j+n_l+1)=0;
        end
    end
end
%-------------------高斯内核与图像卷积-------------------------------------
g_LoG=zeros(M,N);
%对原图进行扩展，方便处理边界
I_pad=padarray(I,[n_l+1,n_l+1],'symmetric');
for i=1:M
    for j=1:N
        %获得图像子块区域
        Block=I_pad(i:i+2*n_l,j:j+2*n_l);
        %用Kirsch内核对子区域卷积     
        g_LoG(i,j)=sum(sum(Block.*Gaussian));
    end
end
n=25;
sigma=4;
type='symmetric';
g_LoG=GaussianBlur(n,I,sigma,type);
% imshow(g_LoG)
%-------------------拉普拉斯算子与上述卷积结果进行卷积---------------------
%拉普拉斯算子
L=[1 1 1;
    1 -8 1;
    1 1 1];
%系数
a=-1;
len=1;
%对原图进行扩展，方便处理边界
g_LoG_pad=padarray(g_LoG,[len+1,len+1],'symmetric');
g_Laplacian=zeros(M,N);
g_L=zeros(M,N);
for i=1+len:M-len
    for j=1+len:N-len
        %从扩展图像中，取出局部图像
        Block=g_LoG_pad(i-len:i+len,j-len:j+len);
        %保留拉普拉斯算子的运算结果
        g_L(i-len,j-len)=a*sum(sum(Block.*L));
    end
end
%计算Laplacian算子标定值
m_min=min(g_L(:));
img=(-m_min+g_L);
g_BiaoDing=img/max(img(:));
% imshow(g_BiaoDing)
 
%-----------------------0值跨界边界搜索------------------------------------
%设定门限
% th=0.04;
th=0.112;
n_max=max(g_L(:));
T=th*n_max;
ind=find(g_L>T)
%图像二值化
g_b=zeros(M,N);
g_b(ind)=1;
 
%对二值图像周边以0扩展
g_b_pad=padarray(g_b,[1,1]);
%0值跨界,寻找边界
g_z=zeros(M,N);
for i=2:M+1
    for j=2:N+1
        %以0或者1为中心
        if g_b_pad(i,j)==1
            %计算以(i,j)为中心点的四个对角是否相异
            if xor(g_b_pad(i-1,j-1),g_b_pad(i+1,j+1))||... %左对角
                xor(g_b_pad(i-1,j+1),g_b_pad(i+1,j-1))||...%右对角
                xor(g_b_pad(i,j-1),g_b_pad(i,j+1))||...%行相对
                xor(g_b_pad(i-1,j),g_b_pad(i+1,j)) %列相对
                g_z(i-1,j-1)=1;
            end
        end
    end
end
imshow(g_z)

高斯模糊GaussianBlur函数：

 
function [g]=GaussianBlur(n,I,sigma,type)
 
    [M,N]=size(I);
 
    %生成高斯核函数
 
    G=GaussianKernelG(n,sigma);
 
    %平滑图像
 
    n_l=n-1;
 
    g=zeros(M,N);
 
    %对原图进行扩展，方便处理边界
 
    I_pad=padarray(I,[n_l,n_l],type);
 
    for i=1:M
 
        for j=1:N
 
            %获得图像子块区域
 
            Block=I_pad(i:i+n_l,j:j+n_l);
 
            %用Kirsch内核对子区域卷积     
 
            g(i,j)=sum(sum(Block.*G));
 
        end
 
    end
 
    %归一化
 
    g=g/max(g(:));
 
end
 
%生成高斯核函数
 
% n 核函数的大小
 
function G=GaussianKernelG(n,sigma)
 
    n_l=floor(n/2);
 
    %初始化
 
    G=zeros(n,n);
 
    %产生高斯核矩阵
 
    for i=-n_l:n_l
 
        for j=-n_l:n_l  
 
            d=i^2+j^2;
 
            G(i+n_l+1,j+n_l+1)=exp(-(d)/(2*sigma^2));
 
        end
 
    end
 
    %寻找最小值
 
    m=sum(G(:));
 
    %取整
 
    G=G/m;
 
    %将大于3*delta的取值置零
 
    for i=-n_l:n_l
 
        for j=-n_l:n_l 
 
            d=sqrt(i^2+j^2);
 
            if d>3*sigma
 
                G(i+n_l+1,j+n_l+1)=0;
 
            end
 
        end
 
    end
 
end