01年数学建模A题程序

程序1：计算最大内切圆半径的粗略范围。运行前请先将程序中“F:\01A\图片\0.bmp”改为图片在你电脑中的实际储存路径。运行时用鼠标在切片图中点出内切圆的大致位置；运行后显示的r是半径估计区间的下限，r1是上限。clear,clcima=imread('F:\01A\图片\0.bmp');%读入一张图片figureimshow(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));endwhiles==0r=r+1;[xxyy]=yuan(j,i,r);s=0;fork=1:length(xx)s=s+ima(yy(k),xx(k));endendr=r-1%rminr1=sqrt((sum(sum(ones(size(ima))))-sum(sum(ima)))/pi)程序2：计算最大内切圆半径、圆心坐标，作投影图与三维图。运行前同样要修改路径。运行时要调用函数yuan.m，在命令窗口当tt=100时程序运行完毕，运行时间一般不超过2分钟。运行结束后，工作空间中，aax，aay为利用计算出的数据拟合后的球心横、纵坐标。ar为从每张图计算出的最大内切圆半径。clear,clcr=28;tt=0;forii=0:99tt=tt+1ima=imread(['F:\01A\图片\'int2str(tt-1)'.bmp']);ima1=ima;ima1([1512],:)=0;ima1(:,[1512])=0;fori=3:3:510forj=3:3:510ifsum(sum(ima(i-1:i+1,[j-2j+2])))+sum(sum(ima([i-2i+2],j-1:j+1)))==12ima1(i-1:i+1,j-1:j+1)=0;endendend%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)==1fl=floor(r/sqrt(2));ima(yeg(k)-fl:yeg(k)+fl,xeg(k)-fl:xeg(k)+fl)=1;endfork1=1:length(xx)ima(yy(k1),xx(k1))=1;%剔除区域C以外的点endend[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));iftlen2len2=t;endendiflen2len22len22=len2;z=k;endendlen2=500000;len22=0;fork=-0.5:0.1:0.5fork0=-0.5:0.1:0.5len2=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));iftlen2len2=t;endendiflen2len22len22=len2;kx=k0;ky=k;endendendax(tt)=yc(z)+ky-256;ay(tt)=xc(z)+kx-256;ar(tt)=sqrt(len22);endz=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)gridonfigure(2)plot(aax,z)gridonfigure(3)plot(aay,z)gridonfigure(4)plot3(aax,aay,z,'r')gridon函数yuan.mfunction[xx,yy]=yuan(x,y,r)k=1;forj=x-r+1:x-1fori=y-r+1:y-1if(x-j)*(x-j)+(y-i)*(y-i)r*rxx(k)=j;yy(k)=i;k=k+1;break;endendendfori=y-r+1:y-1forj=x-r+1:x-1if(x-j)*(x-j)+(y-i)*(y-i)r*rxx(k)=j;yy(k)=i;k=k+1;break;endendendxx(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];以上程序根据国防科技大学向为、王瑛、伍微的论文而编写，但在确定最大圆半径时并不是用二分法，请大家注意。