01年数学建模A题程序.docVIP

程序1：计算最大内切圆半径的粗略范围。 运行前请先将程序中“F:\01A\图片\0.bmp”改为图片在你电脑中的实际储存路径。 运行时用鼠标在切片图中点出内切圆的大致位置； 运行后显示的r是半径估计区间的下限，r1是上限。 clear,clc ima=imread(F:\01A\图片\0.bmp); %读入一张图片 figure imshow(ima) %显示这张图像 [j,i]=ginput(1); %在坐标上取圆心 i=round(i);j=round(j); r=1;s=0; [xx yy]=yuan(j,i,r); for k=1:length(xx) s=s+ima(yy(k),xx(k)); end while s==0 r=r+1; [xx yy]=yuan(j,i,r); s=0; for k=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：计算最大内切圆半径、圆心坐标，作投影图与三维图。 运行前同样要修改路径。 运行时要调用函数yuan.m，在命令窗口当tt=100时程序运行完毕，运行时间一般不超过2分钟。 运行结束后，工作空间中，aax，aay为利用计算出的数据拟合后的球心横、纵坐标。ar为从每张图计算出的最大内切圆半径。 clear,clc r=28;tt=0; for ii=0:99 tt=tt+1 ima=imread([F:\01A\图片\ int2str(tt-1) .bmp]); ima1=ima; ima1([1 512],:)=0;ima1(:,[1 512])=0; for i=3:3:510 for j=3:3:510 if sum(sum(ima(i-1:i+1,[j-2 j+2])))+sum(sum(ima([i-2 i+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); for k=1:length(xeg) [xx yy]=yuan(xeg(k),yeg(k),r); if rem(k,5)==1 fl=floor(r/sqrt(2)); ima(yeg(k)-fl:yeg(k)+fl,xeg(k)-fl:xeg(k)+fl)=1; end for k1=1:length(xx) ima(yy(k1),xx(k1))=1; %剔除区域C以外的点 end end [yc,xc]=find(1-ima); len2=500000;len22=0; for k=1:length(xc) len2=500000; for k1=1:length(xeg) t=(xc(k)-xeg(k1))*(xc(k)-xeg(k1))+(yc(k)-yeg(k1))*(yc(k)-yeg(k1)); if tlen2 len2=t; end end if len2len22 len22=len2; z=k; end end len2=500000;len22=0; for k=-0.5:0.1:0.5 for k0=-0.5:0.1:0.5 len2=500000; for k1=1:length(xeg) xz=xc(z)+k0;yz=yc(z)+k; t=(xz-xeg(k1))*(xz-xeg(k1))+(yz-yeg(k1))*(yz-yeg(k1)); if tlen2 len2=t; end end if len2len22