图像分割——边缘检测——LoG（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