【图像处理】基于灰度矩的亚像素边缘检测方法理论及MATLAB实现

1 基于灰度矩的亚像素边缘检测理论

参考文献：亚像素边缘检测技术研究_张美静

2 MATLAB实现

2.1 main.m

clear;
clc;
tic;%启动计时器，计算程序运行时间
tau=25;
delt=0.5;
N=7;
picture_init=imread('Pic1_2.bmp');
subplot(221);
imshow(picture_init),title('原图像'); 
picture_double=double(picture_init);
[height,wideth]=size(picture_double);

index=1;
fid=fopen('SubpixelEdgeData.txt', 'w');
yaxiangsu_e=zeros(height,wideth);%定义矩阵，并初始化为0，用于边缘的图像的显示
for j=4:1:height-3
    for i=4:1:wideth-3
         m1=conv(picture_double,1,j,i);%计算一阶灰度距
         m2=conv(picture_double,2,j,i);%计算二阶灰度距
         sigma=sqrt(m2-m1^2);
         if sigma>tau
             m3=conv(picture_double,3,j,i);%计算三阶灰度距
             s=(m3+2*m1^3-m1*m2*3)/(sigma^3);
             p1=(1+s*sqrt(1.0/(4+s^2)))/2;%归一化参数
             p2=1-p1;
             h1=m1-sigma*sqrt(p2/p1);
             h2=m1+sigma*sqrt(p1/p2);
             if abs(h1-h2)>sigma*2
                 A=min(p1,p2);
                 x=fzero(@(x)x-0.5*sin(2*x)-A*pi,1.42);%求解超越方程，得到x的值
                 rou =cos(x);
                 if rou <=delt*2/N
                     [x0,y0]=zhongxin(picture_double,j,i);%计算模板圆的灰度重心坐标
                     sin_o=y0/sqrt(x0^2+y0^2);
                     cos_o=x0/sqrt(x0^2+y0^2);
                     Xs(index)=i+rou*cos_o*N/2;%计算亚像素横坐标x坐标值,保存在数组Xs中，方便查看
                     Ys(index)=j+rou*sin_o*N/2;%计算亚像素纵坐标y坐标值
                     fprintf(fid, '%f\t%f\r\n', Xs(index), Ys(index));%将得到亚像素边缘数据保存到txt文件中方便查看
                     Xs_int=round(Xs(index));%取整用于显示结果
                     Ys_int=round(Ys(index));
                     yaxiangsu_e(Ys_int,Xs_int)=1;
                     index=index+1;
                 end
             end
         end      
    end
end
subplot(222);
imshow(yaxiangsu_e),title('灰度矩亚像素边缘检测结果'); 

subplot(223);
I41=imfill(yaxiangsu_e,'holes');
imshow(I41)
title('孔洞填充图像');
% 提取最外围边缘
subplot(224);
I4=bwperim(I41);                   
imshow(I4); title('边缘图像');
% 去除面积小于150px物体
% subplot(224); 
% I5=bwareaopen(I4,100);    
% imshow(I5);

2.2 Conv.m

function result=conv(picture,nsqure,j,i)
%-----------------------计算模板与灰度矩阵的卷积---------------------%
M=[0 0.00913767235 0.021840193 0.025674188 0.021840193 0.00913767235 0;0.00913767235 0.025951560 0.025984481 0.025984481 0.025984481 0.025951560   0.00913767235 ;0.021840193   0.025984481   0.025984481 0.025984481 0.025984481 0.025984481   0.021840193; 0.025674188   0.025984481   0.025984481 0.025984481 0.025984481 0.025984481   0.025674188;    0.021840193   0.025984481   0.025984481 0.025984481 0.025984481 0.025984481   0.021840193;    0.00913767235 0.025951560   0.025984481 0.025984481 0.025984481 0.025951560   0.00913767235; 0 0.00913767235 0.021840193 0.025674188 0.021840193 0.00913767235 0];%卷积模板

result=picture(j-3,i-3)^nsqure*M(1)+picture(j-2,i-3)^nsqure*M(2)+picture(j-1,i-3)^nsqure*M(3)+picture(j,i-3)^nsqure*M(4)+picture(j+1,i-3)^nsqure*M(5)+picture(j+2,i-3)^nsqure*M(6)+picture(j+3,i-3)^nsqure*M(7)+picture(j-3,i-2)^nsqure*M(8)+picture(j-2,i-2)^nsqure*M(9)+picture(j-1,i-2)^nsqure*M(10)+picture(j,i-2)^nsqure*M(11)+picture(j+1,i-2)^nsqure*M(12)+picture(j+2,i-2)^nsqure*M(13)+picture(j+3,i-2)^nsqure*M(14)+picture(j-3,i-1)^nsqure*M(15)+picture(j-2,i-1)^nsqure*M(16)+picture(j-1,i-1)^nsqure*M(17)+picture(j,i-1)^nsqure*M(18)+picture(j+1,i-1)^nsqure*M(19)+picture(j+2,i-1)^nsqure*M(20)+picture(j+3,i-1)^nsqure*M(21)+picture(j-3,i)^nsqure*M(22)+picture(j-2,i)^nsqure*M(23)+picture(j-1,i)^nsqure*M(24)+picture(j,i)^nsqure*M(25)+picture(j+1,i)^nsqure*M(26)+picture(j+2,i)^nsqure*M(27)+picture(j+3,i)^nsqure*M(28)+picture(j-3,i+1)^nsqure*M(29)+picture(j-2,i+1)^nsqure*M(30)+picture(j-1,i+1)^nsqure*M(31)+picture(j,i+1)^nsqure*M(32)+picture(j+1,i+1)^nsqure*M(33)+picture(j+2,i+1)^nsqure*M(34)+picture(j+3,i+1)^nsqure*M(35)+picture(j-3,i+2)^nsqure*M(36)+picture(j-2,i+2)^nsqure*M(37)+picture(j-1,i+2)^nsqure*M(38)+picture(j,i+2)^nsqure*M(39)+picture(j+1,i+2)^nsqure*M(40)+picture(j+2,i+2)^nsqure*M(41)+picture(j+3,i+2)^nsqure*M(42)+picture(j-3,i+3)^nsqure*M(43)+picture(j-2,i+3)^nsqure*M(44)+picture(j-1,i+3)^nsqure*M(45)+picture(j,i+3)^nsqure*M(46)+picture(j+1,i+3)^nsqure*M(47)+picture(j+2,i+3)^nsqure*M(48)+picture(j+3,i+3)^nsqure*M(49);

2.3 Zhongxin.m

function [x,y]=zhongxin(picture,j,i)
%-------用于计算灰度圆的灰度重心坐标-------%
huidu_sum=0;
huidu_x=0;
huidu_y=0;
for m=j-3:1:j+3
    for n=i-3:1:i+3
        huidu_sum=huidu_sum+picture(m,n);
        huidu_x=huidu_x+picture(m,n)*n;
        huidu_y=huidu_y+picture(m,n)*m;
    end
end
huidu_sum=huidu_sum-picture(j-3,i-3)-picture(j-3,i+3)-picture(j+3,i-3)-picture(j+3,i+3);
huidu_x=huidu_x-picture(j-3,i-3)*(i-3)-picture(j-3,i+3)*(i+3)-picture(j+3,i-3)*(i-3)-picture(j+3,i+3)*(i+3);
huidu_y=huidu_y-picture(j-3,i-3)*(j-3)-picture(j-3,i+3)*(j-3)-picture(j+3,i-3)*(j+3)-picture(j+3,i+3)*(j+3);
x=huidu_x/(huidu_sum+eps);
y=huidu_y/(huidu_sum+eps);

实验结果图