01年数学建模A题程序文档格式.docx

01年数学建模A题程序文档格式.docx

《01年数学建模A题程序文档格式.docx》由会员分享，可在线阅读，更多相关《01年数学建模A题程序文档格式.docx（6页珍藏版）》请在冰豆网上搜索。

imshow（ima）%显示这张图像

[j,i]=ginput

（1）;

%在坐标上取圆心

i=round（i）;

j=round（j）;

r=1;

s=0;

[xxyy]=yuan（j,i,r）;

fork=1:

length（xx）

s=s+ima（yy（k）,xx（k））;

end

whiles==0

r=r+1;

[xxyy]=yuan（j,i,r）;

s=0;

fork=1:

end

r=r-1%rmin

r1=sqrt（（sum（sum（ones（size（ima））））-sum（sum（ima）））/pi）

程序2：

计算最大内切圆半径、圆心坐标，作投影图与三维图。

运行前同样要修改路径。

运行时要调用函数，在命令窗口当tt=100时程序运行完毕，运行时间一般不超过2分钟。

运行结束后，工作空间中，aax，aay为利用计算出的数据拟合后的球心横、纵坐标。

ar为从每张图计算出的最大内切圆半径。

r=28;

tt=0;

forii=0:

99

tt=tt+1

ima=imread（['

int2str（tt-1）'

.bmp'

]）;

ima1=ima;

ima1（[1512],:

）=0;

ima1（:

[1512]）=0;

fori=3:

3:

510

forj=3:

ifsum（sum（ima（i-1:

i+1,[j-2j+2]）））+sum（sum（ima（[i-2i+2],j-1:

j+1）））==12

ima1（i-1:

i+1,j-1:

j+1）=0;

%eg=edge（ima,'

sobel'

求出切片的边缘点的位置

%[yeg,xeg]=find（eg）;

边缘点的横纵坐标

[yeg,xeg]=find（ima1）;

length（xeg）

[xxyy]=yuan（xeg（k）,yeg（k）,r）;

ifrem（k,5）==1

fl=floor（r/sqrt

（2））;

ima（yeg（k）-fl:

yeg（k）+fl,xeg（k）-fl:

xeg（k）+fl）=1;

fork1=1:

ima（yy（k1）,xx（k1））=1;

%剔除区域C以外的点

[yc,xc]=find（1-ima）;

len2=500000;

len22=0;

length（xc）

len2=500000;

t=（xc（k）-xeg（k1））*（xc（k）-xeg（k1））+（yc（k）-yeg（k1））*（yc（k）-yeg（k1））;

ift<

len2

len2=t;

iflen2>

len22

len22=len2;

z=k;

fork=:

:

fork0=:

xz=xc（z）+k0;

yz=yc（z）+k;

t=（xz-xeg（k1））*（xz-xeg（k1））+（yz-yeg（k1））*（yz-yeg（k1））;

kx=k0;

ky=k;

ax（tt）=yc（z）+ky-256;

ay（tt）=xc（z）+kx-256;

ar（tt）=sqrt（len22）;

z=0:

99;

A=polyfit（z,ax,7）;

B=polyfit（z,ay,7）;

aax=polyval（A,z）;

aay=polyval（B,z）;

figure

（1）

plot（aax,aay）

gridon

figure

（2）

plot（aax,z）

figure（3）

plot（aay,z）

figure（4）

plot3（aax,aay,z,'

r'

）

函数

function[xx,yy]=yuan（x,y,r）

k=1;

forj=x-r+1:

x-1

fori=y-r+1:

y-1

if（x-j）*（x-j）+（y-i）*（y-i）<

r*r

xx（k）=j;

yy（k）=i;

k=k+1;

break;

fori=y-r+1:

forj=x-r+1:

xx（k:

k+1）=[x-r+1x];

yy（k:

k+1）=[yy-r+1];

xx=[xxxx2*x-xx2*x-xx];

yy=[yy2*y-yyyy2*y-yy];

以上程序根据国防科技大学向为、王瑛、伍微的论文而编写，但在确定最大圆半径时并不是用二分法，请大家注意。