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);
实验结果图
