Zernike矩亚像素边缘检测

function zernike7(I) I=imread('2.jpg'); % 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] imshow(J) % hold on %显示针尖的位置 % plot(y1sub,1:m,'r') % plot(y2sub,1:m,'r') % ymin=min(Z(find(Z(:,2)>696),2)) % Z(find(Z(:,2)==ymin),:) % % ymax=max(Z(find(Z(:,2)<696),2)) % Z(find(Z(:,2)==ymax),:) % % imshow(J) % hold on % plot(Z(find(Z(:,2)==ymin),2),1:m,'r') % plot(Z(find(Z(:,2)==ymax),2),1:m,'r')