matlab程序.txt

2008年度全国数学建模竞赛用程序

%得二值图，边缘检测，获轮廓%
clear all
close all
 BW=imread('a.jpg');
 I=im2bw(BW);
 I= bwperim(I,8);   % 直接给定边界轮廓,通过八方位链码跟踪%
 I(1,:)=0;
 I(768,:)=0;
 I(:,1)=0;
 I(:,1024)=0;
 for i=120:260
for j=1:375
I1(i,j)=I(i,j);
end
end
 for i=120:260
  for j=375:500
I2(i,j)=I(i,j);
end
end
 for i=120:260
  for j=540:700
  I3(i,j)=I(i,j);
end
end
 for i=430:550
 for j=200:375
  I4(i,j)=I(i,j);
end
end
for i=430:550
for j=530:650
I5(i,j)=I(i,j);
end
end


n=0;
for i=120:260
for j=1:375
if I1(i,j)==1
n=n+1;
end
end
end
c1=([linspace(0,0,n);linspace(0,0,n)]);
k=1;
for i=120:260
for j=1:375
if I1(i,j)==1,
c1(1,k)=i;c1(2,k)=j;k=k+1;
end
end
end


n=0;
 for i=120:260
  for j=375:500
if I2(i,j)==1
n=n+1;
end
end
end
c2=([linspace(0,0,n);linspace(0,0,n)]);
k=1;
 for i=120:260
  for j=375:500
if I2(i,j)==1,
c2(1,k)=i;c2(2,k)=j;k=k+1;
end
end
end

n=0;
for i=120:260
  for j=540:700
if I3(i,j)==1
n=n+1;
end
end
end
c3=([linspace(0,0,n);linspace(0,0,n)]);
k=1;
for i=120:260
  for j=540:700
if I3(i,j)==1,
c3(1,k)=i;c3(2,k)=j;k=k+1;
end
end
end


n=0;
 for i=430:550
 for j=200:375
if I4(i,j)==1
n=n+1;
end
end
end
c4=([linspace(0,0,n);linspace(0,0,n)]);
k=1;
 for i=430:550
 for j=200:375
if I4(i,j)==1,
c4(1,k)=i;c4(2,k)=j;k=k+1;
end
end
end

n=0;
for i=430:550
for j=530:650
if I5(i,j)==1
n=n+1;
end
end
end
c5=([linspace(0,0,n);linspace(0,0,n)]);
k=1;
for i=430:550
for j=530:650
if I5(i,j)==1,
c5(1,k)=i;c5(2,k)=j;k=k+1;
end
end
end

%切线法求坐标%
m=min(c1(1,:));k=max(c1(1,:));a1=(m+k)/2;
n=min(c1(2,:));l=max(c1(2,:));b1=(n+l)/2;
m=min(c2(1,:));k=max(c2(1,:));a2=(m+k)/2;
n=min(c2(2,:));l=max(c2(2,:));b2=(n+l)/2;
m=min(c3(1,:));k=max(c3(1,:));a3=(m+k)/2;;
n=min(c3(2,:));l=max(c3(2,:));b3=(n+l)/2;
m=min(c4(1,:));k=max(c4(1,:));a4=(m+k)/2;
n=min(c4(2,:));l=max(c4(2,:));b4=(n+l)/2;
m=min(c5(1,:));k=max(c5(1,:));a5=(m+k)/2;
n=min(c5(2,:));l=max(c5(2,:));b5=(n+l)/2;
S(1,1)=a1;S(1,2)=a2;S(1,3)=a3;S(1,4)=a4;S(1,5)=a5;
S(2,1)=b1;S(2,2)=b2;S(2,3)=b3;S(2,4)=b4;S(2,5)=b5;S=S';

%利用椭圆拟合法求坐标%
[A,B,Xc,Yc,P]=ellipsefit(c1(1,:),c1(2,:));
R(1,1)=Xc;R(1,2)=Yc;
[A,B,Xc,Yc,P]=ellipsefit(c2(1,:),c2(2,:));
R(2,1)=Xc;R(2,2)=Yc;
[A,B,Xc,Yc,P]=ellipsefit(c3(1,:),c3(2,:));
R(3,1)=Xc;R(3,2)=Yc;
[A,B,Xc,Yc,P]=ellipsefit(c4(1,:),c4(2,:));
R(4,1)=Xc;R(4,2)=Yc;
[A,B,Xc,Yc,P]=ellipsefit(c5(1,:),c5(2,:));
R(5,1)=Xc;R(5,2)=Yc;
C2(:,1)=R(:,2);C2(:,2)=-R(:,1);
T=[-512,384;-512,384;-512,384;-512,384;-512,384];
C2=C2+T;C2=C2/3.78,R
C1(:,1)=S(:,2);C1(:,2)=-S(:,1);
T=[-512,384;-512,384;-512,384;-512,384;-512,384];
C1=C1+T;C1=C1/3.78,S


A=C';
h=417.2;
for i=1:5;for j=1:5;cos(i,j)=(h^2+A(1,i)*A(1,j)+A(2,i)*A(2,j))/sqrt((h^2+A(1,i)^2+A(2,i)^2)*(h^2+A(1,j)^2+A(2,j)^2));end;end,cos

[x,y,z]=solve('a3*a5*x+(y*a2+z*a1)*a4*a5-(y*a1-z*a2)*a6-a8=0','x*a3*a6+(y*a2+z*a1)*a4*a6+(y*a1-z*a2)*a5-a9=0','-x*a4+(y*a2+z*a1)*a3-a7=0')
 
a=[ 0.9264086,0.3765197,0.9378301,0.3470948,0.9994888,-0.0319720,487.2626,15.00192,-8.177267];
a1=a(1);a2=a(2);a3=a(3);a4=a(4);a5=a(5);a6=a(6);a7=a(7);a8=a(8);a9=a(9);