matlab最大内切圆程序---图像处理

I=imread('D:/xx.bmp');threshold=graythresh(I);%计算将灰度图像转化为二值图像所需的门限I1=im2bw(I,threshold);%将灰度图像转化为二值图像I2=~I1;dim=size(I1);col=round(dim(2)/2)-90;%计算起始点列坐标row=find(I1(:,col),1);%计算起始点行坐标?connectivity=8;num_points=1390;contour=bwtraceboundary(I1,[row,col],'N',connectivity,num_points);%提取边界b=edge(I2,'sobel');%findingedges[height,width]=size(b);R=60;fori=60:1:width-60forj=60:1:height-60forr=60:1:120fora=0:1:360t=1;x=i+r*sin(a);y=j+r*cos(a);x=round(x);y=round(y);if(x0)&&(xwidth)&&(y0)&&(yheight)ifb(y,x)==1||b(j,i)==1t=0;break;endendendift==0break;endifRrR=r;ox=i;oy=j;endendendendi=1;forc=pi/10000:pi/10000:2*pi;x(i)=ceil(cos(c)*R+ox);y(i)=ceil(sin(c)*R+oy);i=i+1;endpx=I2;py=I;ford=1:(i-3)if(x(d)500)%不超出范围if(y(d)500)if(x(d)0)if(y(d)0)px(y(d),x(d))=1;%循环画圆py(y(d),x(d))=1;endendendendendfigure,imshow(px);title('yuan1');figure,imshow(py);title('yuan2');