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

1 边缘检测理论

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

2 MATLAB实现

function zernike7(I)

I=imread('Pic1_3.bmp');
subplot(221)
imshow(I)
title('原图像');
% 7*7Zernike模板
M00=...
[
       0    0.0287    0.0686    0.0807    0.0686    0.0287    0
    0.0287    0.0815    0.0816    0.0816    0.0816    0.0815    0.0287
    0.0686    0.0816    0.0816    0.0816    0.0816    0.0816    0.0686
    0.0807    0.0816    0.0816    0.0816    0.0816    0.0816    0.0807
    0.0686    0.0816    0.0816    0.0816    0.0816    0.0816    0.0686
    0.0287    0.0815    0.0816    0.0816    0.0816    0.0815    0.0287
       0    0.0287    0.0686    0.0807    0.0686    0.0287    0
];
M11R=...
[
       0    -0.015    -0.019    0    0.019    0.015      0
    -0.0224    -0.0466    -0.0233    0    0.0233    0.0466    0.0224
    -0.0573    -0.0466    -0.0233    0    0.0233    0.0466    0.0573
    -0.069    -0.0466    -0.0233    0    0.0233    0.0466    0.069
    -0.0573    -0.0466    -0.0233    0    0.0233    0.0466    0.0573
    -0.0224    -0.0466    -0.0233    0    0.0233    0.0466    0.0224
       0    -0.015    -0.019    0    0.019    0.015      0
];
M11I=...
[
       0    -0.0224    -0.0573    -0.069    -0.0573    -0.0224      0
    -0.015    -0.0466    -0.0466    -0.0466    -0.0466    -0.0466    -0.015
    -0.019    -0.0233    -0.0233    -0.0233    -0.0233    -0.0233    -0.019
       0       0       0       0       0       0      0
    0.019    0.0233    0.0233    0.0233    0.0233    0.0233    0.019
    0.015    0.0466    0.0466    0.0466    0.0466    0.0466    0.015
       0    0.0224    0.0573    0.069    0.0573    0.0224      0
];
M20=...
[
       0    0.0225    0.0394    0.0396    0.0394    0.0225    0
    0.0225    0.0271    -0.0128    -0.0261    -0.0128    0.0271    0.0225
    0.0394    -0.0128    -0.0528    -0.0661    -0.0528    -0.0128    0.0394
    0.0396    -0.0261    -0.0661    -0.0794    -0.0661    -0.0261    0.0396
    0.0394    -0.0128    -0.0528    -0.0661    -0.0528    -0.0128    0.0394
    0.0225    0.0271    -0.0128    -0.0261    -0.0128    0.0271    0.0225
       0    0.0225    0.0394    0.0396    0.0394    0.0225    0
];


if length(size(I))==3 I=rgb2gray(I); end
I=im2bw(I,0.6);
K=double(I);
[m n]=size(K);
xs=double(zeros(m,n));
ys=double(zeros(m,n));

% 卷积运算
A11I=conv2(M11I,K);
A11R=conv2(M11R,K);
A20=conv2(M20,K);

% 截掉多余部分
A11I=A11I(4:end-3,4:end-3);
A11R=A11R(4:end-3,4:end-3);
A20=A20(4:end-3,4:end-3);

J=zeros(size(K));
boundary=J;
theta=atan2(A11I,A11R);%计算theta

%计算边缘的三个参数
A11C=A11R.*cos(theta)+A11I.*sin(theta);
l=A20./A11C;
k=1.5*A11C./((1-l.^2).^1.5);
e=abs(l)>1/3.5;
k(e)=0;

%边缘判断条件
a=abs(l)<1/sqrt(2)*2/7;
b=abs(k)>max(I(:))/10;
% a,b分别为距离和边缘强度判断结果
J(a&b)=1;
%将图像的最边缘去除
% boundary(2:end-1,2:end-1)=1;
% J(~boundary)=0;
format short
% [x,y]=find(J==1);%边缘的像素级坐标
% O=[x y];
% Z=[x+l(find(J==1)).*cos(theta(find(J==1))) y+l(find(J==1)).*sin(theta(find(J==1)))];%亚像素坐标
% % fprintf('%.4f %.4f\n',Z');

[L,num]=bwlabel(J,8);%对二值图像进行标记
%自动化搜索连通域
s=zeros(1,num);
for i=1:num
   s(i)=size(find(L==i),1);
end
[bwL,label]=sort(s,'descend');
if label(1)<label(2)
    index1=label(1);
    index2=label(2);
else
    index1=label(2);
    index2=label(1);
end

%计算左边探针的最前端坐标
[r1,c1]=find(L==index1);
A1=[r1 c1];
y1=max(A1(:,2));%该连通域中y最大值为针尖处
x1=max(A1(find(A1(:,2)==y1),1));
x1sub=x1+3.5*l(x1,y1)*cos(theta(x1,y1));
y1sub=y1+3.5*l(x1,y1)*sin(theta(x1,y1));

%计算最右边探针的最前端坐标
[r2,c2]=find(L==index2);
A2=[r2 c2];
y2=min(A2(:,2));%该连通域中y最小为连通域
x2=max(A2(find(A2(:,2)==y2),1));
x2sub=x2+3.5*l(x2,y2)*cos(theta(x2,y2));
y2sub=y2+3.5*l(x2,y2)*sin(theta(x2,y2));

% [x1sub y1sub],[x2sub,y2sub]
subplot(222)
imshow(J)
title('检测结果图像');
% figure;
% imcontour(J,1)

% 边界提取
% subplot(122);
% bwimg = bwmorph(J,'remove');
% imshow(bwimg)

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

实验图