01年数学建模A题程序.docx
《01年数学建模A题程序.docx》由会员分享，可在线阅读，更多相关《01年数学建模A题程序.docx（6页珍藏版）》请在冰豆网上搜索。
01年数学建模A题程序
程序1：
计算最大内切圆半径的粗略范围。
运行前请先将程序中“F:
\01A\图片\”改为图片在你电脑中的实际储存路径。
运行时用鼠标在切片图中点出内切圆的大致位置；
运行后显示的r是半径估计区间的下限，r1是上限。
clear,clc
ima=imread（'F:
\01A\图片\'）;%读入一张图片
figure
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:
length（xx）
s=s+ima（yy（k）,xx（k））;
end
end
r=r-1%rmin
r1=sqrt（（sum（sum（ones（size（ima））））-sum（sum（ima）））/pi）
程序2：
计算最大内切圆半径、圆心坐标，作投影图与三维图。
运行前同样要修改路径。
运行时要调用函数，在命令窗口当tt=100时程序运行完毕，运行时间一般不超过2分钟。
运行结束后，工作空间中，aax，aay为利用计算出的数据拟合后的球心横、纵坐标。
ar为从每张图计算出的最大内切圆半径。
clear,clc
r=28;tt=0;
forii=0:
99
tt=tt+1
ima=imread（['F:
\01A\图片\'int2str（tt-1）'.bmp']）;
ima1=ima;
ima1（[1512],:
）=0;ima1（:
[1512]）=0;
fori=3:
3:
510
forj=3:
3:
510
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;
end
end
end
%eg=edge（ima,'sobel'）;求出切片的边缘点的位置
%[yeg,xeg]=find（eg）;边缘点的横纵坐标
[yeg,xeg]=find（ima1）;
fork=1:
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;
end
fork1=1:
length（xx）
ima（yy（k1）,xx（k1））=1;%剔除区域C以外的点
end
end
[yc,xc]=find（1-ima）;
len2=500000;len22=0;
fork=1:
length（xc）
len2=500000;
fork1=1:
length（xeg）
t=（xc（k）-xeg（k1））*（xc（k）-xeg（k1））+（yc（k）-yeg（k1））*（yc（k）-yeg（k1））;
iftlen2=t;
end
end
iflen2>len22
len22=len2;
z=k;
end
end
len2=500000;len22=0;
fork=:
:
fork0=:
:
len2=500000;
fork1=1:
length（xeg）
xz=xc（z）+k0;yz=yc（z）+k;
t=（xz-xeg（k1））*（xz-xeg（k1））+（yz-yeg（k1））*（yz-yeg（k1））;
iftlen2=t;
end
end
iflen2>len22
len22=len2;kx=k0;ky=k;
end
end
end
ax（tt）=yc（z）+ky-256;
ay（tt）=xc（z）+kx-256;
ar（tt）=sqrt（len22）;
end
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）
gridon
figure（3）
plot（aay,z）
gridon
figure（4）
plot3（aax,aay,z,'r'）
gridon
函数
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）xx（k）=j;
yy（k）=i;
k=k+1;
break;
end
end
end
fori=y-r+1:
y-1
forj=x-r+1:
x-1
if（x-j）*（x-j）+（y-i）*（y-i）xx（k）=j;
yy（k）=i;
k=k+1;
break;
end
end
end
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];
以上程序根据国防科技大学向为、王瑛、伍微的论文而编写，但在确定最大圆半径时并不是用二分法，请大家注意。
升级会员 