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