全国大学生数学建模竞赛2001a赛题matlab代码分享

话不多说，直接上代码

for num_graph=0:5:99    %每隔5张图片画一个边界图，这样看得更清楚
    s1='C:\Users\86135\Desktop\competition\2001\2001aphoto\'
    s2=int2str(num_graph)
    s3='.bmp'
    s=strcat(s1,s2,s3)
    graphN=imread(s)    %读入位图
    edgeN=edge(graphN,'log')    %求出边界点
    for i =1:512
        for j=1:512
            if edgeN(i,j)==1
                plot3(j,512-i,num_graph,'.')
                hold on
            end
        end
    end
end
axis equal
box on
tic    %记录运行时间
mat_maxR=[]
mat_meanC0=[]
for l=0:99
    graphN=[]
    edgeN=[]
    x=[]
    y=[]
    co=[]
    s1='C:\Users\86135\Desktop\competition\2001\2001aphoto\'
    s2=int2str(l)
    s3='.bmp'
    s=strcat(s1,s2,s3)
    graphN=imread(s)
    edgeN=edge(graphN,'log')
    [x,y]=find(edgeN==1)
    leftLine=min(y)+28
    topLine=min(x)+28
    rightLine=max(y)-28
    bottomLine=max(x)-28
    maxR=0;
    meanC0=[0;0]
    for r=28:32
        r2=r^2
        co=[]
        isfind02=0    %isfind02是否为最大内切圆
        for i =topLine:bottomLine    %从上到下
            for j=leftLine:rightLine    %从左到右搜索圆心坐标
                isfind0=1
                if graphN(i,j)==0
                    for k=1:length(x)
                        dist=(x(k)-i)^2+(y(k)-j)^2
                        if dist<r2    %当边界点在圆内，则跳出循环
                            isfind0=0    %isfind0为0表示(i,j)不是圆心
                            break
                        end
                    end
                    if isfind0==1    %isfind0为1表示(i,j)是圆心
                        isfind02=1    %设置标识，找到内切圆
                        co=[co,[i;j]]    %存下数据
                    end
                end
            end
        end
        if isfind02==1    %如果找到内切圆，则记下圆心和半径值
            meanC0=round([mean(co(1,:));mean(co(2,:))])
            maxR=r
        else
            break    %如果不是最大内切圆，则停止搜索
        end
    end
    mat_meanC0=[meanC0,mat_meanC0]    %存储全部的圆心坐标
    mat_maxR=[mat_maxR,maxR]    %存储全部的100个半径值
end
toc    %记录运行时间
m=[257;96]
n=[244;445]    %增加两个位点，使三次B样条曲线过初始点和终点
mat_meanC0=[m mat_meanC0 n]
for i =1:99
    x0=mat_meanC0(1,i)    
    x1=mat_meanC0(1,i+1)
    x2=mat_meanC0(1,i+2)
    x3=mat_meanC0(1,i+3)
    y0=mat_meanC0(2,i)
    y1=mat_meanC0(2,i+1)
    y2=mat_meanC0(2,i+2)
    y3=mat_meanC0(2,i+3)
    z0=i-2
    z1=i-1
    z2=i
    z3=i+1
    a0=(x0+4*x1+x2)/6    %设置a0,a1,a2,a3的数值
    a1=(x2-x0)/2
    a2=(x0-2*x1+x2)/2
    a3=(x0-3*x1+3*x2-x3)/(-6)
    b0=(y0+4*y1+y2)/6
    b1=(y2-y0)/2
    b2=(y0-2*y1+y2)/2
    b3=(y0-3*y1+3*y2-y3)/(-6)
    c0=(z0+4*z1+z2)/6
    c1=(z2-z0)/2
    c2=(z0-2*z1+z2)/2
    c3=(z0-3*z1+3*z2-z3)/(-6)
    Bx=[]
    By=[]
    Bz=[]
    for t=0:0.1:1
        bx=a0+a1*t+a2*t^2+a3*t^3
        by=b0+b1*t+b2*t^2+b3*t^3
        bz=c0+c1*t+c2*t^2+c3*t^3
        Bx=[Bx bx]
        By=[By by]
        Bz=[Bz bz]
    end
    plot3(Bx,By,Bz)    %画出图像
    hold on
end
xlabel('x')
ylabel('y')
zlabel('z')
axis equal

本代码为本人手打，运行结果请自行上机演示。