Im=zeros(m,n); for x=5:m-5; for y=5:n-5;
sum1=I(x,y-4)+I(x,y-2)+I(x,y+2)+I(x,y+4);
sum2=I(x-2,y+4)+I(x-1,y+2)+I(x+1,y-2)+I(x+2,y-4); sum3=I(x-2,y+2)+I(x-4,y+4)+I(x+2,y-2)+I(x+4,y-4); sum4=I(x-2,y+1)+I(x-4,y+2)+I(x+2,y-1)+I(x+4,y-2); sum5=I(x-2,y)+I(x-4,y)+I(x+2,y)+I(x+4,y);
sum6=I(x-4,y-2)+I(x-2,y-1)+I(x+2,y+1)+I(x+4,y+2); sum7=I(x-4,y-4)+I(x-2,y-2)+I(x+2,y+2)+I(x+4,y+4); sum8=I(x-2,y-4)+I(x-1,y-2)+I(x+1,y+2)+I(x+2,y+4); sumi=[sum1,sum2,sum3,sum4,sum5,sum6,sum7,sum8]; summax=max(sumi); summin=min(sumi); summ=sum(sumi); b=summ/8;
if (summax+summin+ 4*I(x,y))> (3*summ/8) sumf = summin; else
sumf =summax; end
if sumf > b Im(x,y)=128; else
Im(x,y)=255; end end end for i=1:m
for j =1:n
Icc(i,j)=Icc(i,j)*Im(i,j); end end
for i=1:m for j =1:n
if (Icc(i,j)==128) Icc(i,j)=0; else
Icc(i,j)=1; end; end end
figure,imshow(double(Icc));title('二值化');
%************************************************************************
36
u=Icc;
[m,n]=size(u) %去空洞和毛刺 for x=2:m-1 for y=2:n-1 if u(x,y)==0
if u(x,y-1)+u(x-1,y)+u(x,y+1)+u(x+1,y)>=3 u(x,y)=1; end
else u(x,y)=u(x,y); end end end
figure,imshow(u) %title('去毛刺') for a=2:m-1 for b=2:n-1 if u(a,b)==1 if
abs(u(a,b+1)-u(a-1,b+1))+abs(u(a-1,b+1)-u(a-1,b))+abs(u(a-1,b)-u(a-1,b-1))+abs(u(a-1,b-1)-u(a,b-1))+abs(u(a,b-1)-u(a+1,b-1))+abs(u(a+1,b-1)-u(a+1,b))+abs(u(a+1,b)-u(a+1,b+1))+abs(u(a+1,b+1)-u(a,b+1))~=1 %去空洞
if(u(a,b+1)+u(a-1,b+1)+u(a-1,b))*(u(a,b-1)+u(a+1,b-1)+u(a+1,b))+(u(a-1,b)+u(a-1,b-1)+u(a,b-1))*(u(a+1,b)+u(a+1,b+1)+u(a,b+1))==0 %去毛刺 u(a,b)=0; end end end end end
figure,imshow(u) %title('去空洞')
%************************************************************************* v=~u;
se=strel('square',3); fo=imopen(v,se);
v=imclose(fo,se); %对图像开操作和闭操作 img=bwmorph(v,'thin',Inf); %对图像进行细化 figure,imshow(img) title('细化图')
附录B 特征点提取代码
function j = P (img, x, y, i) switch (i)
case {1, 9}
j = img(x+1, y);
37
case 2
j = img(x + 1, y-1); case 3
j = img(x, y - 1); case 4
j = img(x - 1, y - 1); case 5
j = img(x - 1, y); case 6
j = img(x - 1, y + 1); case 7
j = img(x, y + 1); case 8
j = img(x + 1, y + 1); end
%point函数
function txy=point(thin) count = 1;
txy(count, :) = [0,0,0];
siz=min(size(thin,1),size(thin,2)); for x=40:siz - 40 for y=40:siz - 40 if (thin(y, x) )
CN = 0; for i = 1:8
CN = CN + abs (P(thin, y, x, i) - P(thin, y, x, i + 1)); end if (CN == 2)
txy(count, :) = [x, y,2]; count = count + 1; end
if (CN == 6)
txy(count, :) = [x, y,6]; count = count + 1; end end end end
for i=1:count - 1 x(i) =txy(i, 1); y(i)= txy(i, 2);
38
end
imshow(double(thin)); hold on; plot(x,y,'.'); %cut函数
function txy=cut(thin,txy) s(8,8)=0; delta(8,8)=0; n=size(txy,1); for i=1:8 for j=1:8
mp{i,j}=thin(1+31*(i-1):31+31*(i-1),1+31*(j-1):31+31*(j-1)); s(i,j)=sum(sum(mp{i,j}))/(31*31); mp{i,j}=(mp{i,j}-s(i,j)).^2; delta(i,j)=sum(sum(mp{i,j})); if delta(i,j)<=70 for k=1:n if
(txy(k,1)>=1+31*(i-1)&&txy(k,1)<=31+31*(i-1)&&txy(k,2)>=1+31*(j-1)&&txy(k,2)<=31+31*(j-1)&&txy(k,3)==2)
txy(k,:)=[0,0,0]; end end
end end end
txy=txy(find(txy(:,1)),:); plot(txy(:,1),txy(:,2),'ro');
附录C 图像特征点代码
%single_point函数
function [pxy2,error]=single_point(txy,r) error=0; x=txy(:,1); y=txy(:,2); n=length(x); d(1:n,1:n)=0;
39
for j=1:n
for i=1:n if (i~=j)
d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2); else
d(i,j)=2*r; end end end
[a,b]=min(d); c=find(a>r);
pxy2=txy(c,:);
pxy2=pxy2(find(pxy2(:,3)==2),:); t=size(pxy2,1); if t==0
error=1 else
plot(x,y,'b.'); hold on
plot(pxy2(:,1),pxy2(:,2),'r.'); end
%walk函数
function [error,a,b]=walk(thin,x0,y0,num) error=0;
thin(y0,x0)=0; t1=0;
for n=1:num if error==1
break; else
x=x0; y=y0;
for x=x0-1:x0+1 if error==1 break; else
for y=y0-1:y0+1
t1=sum(sum(thin(y0-1:y0+1,x0-1:x0+1))); if (t1==0||t1>=2)
40