fractalsurface2.m

分形计算的几个程序原码

%d:\sun\fp 为存放数据文件的路径名;
%s 为存放原始数据文件名变量;
%sp 为存放插值数据文件名变量;
%af 为公式(5)中决定分形插值曲面分形维数(粗糙程度)的自由参数 αn,m的值;
%n, m 分别为 x, y 方向插值结点数 N, M;
%nn, mm 分别为 x, y 方向插值后的结点数;
%x(n×1), y(m×1)分别存放 x, y 方向步长值(Δx), (Δy);
%z(n×m)存放插值结点上的原始数据;
%zz(nn×mm)存放插值后结点上的值;
%a(n×1), b(n×1), c(m×1), d(m×1)分别存放公式(4)中的 an, bn, cm, dm;
%cc(n×m), bb(n×m), dd(n×m), kk(n×m)分别存放公式(6)中的 en,m, fn,m, gn,m, kn,m
% This is a F3PT1.m program (fractal surface interpolation).
bds=input('Input your file name please ','s');
af=input('Input argument af == ');
s=['d:\sun\fp\' bds];
sp=['d:\sun\fp\' bds 'p'];
eval (['load ' s]);
eval(['z=' bds]);
eval(['clear ' bds]);
[m,n]=size(z);
x=0:100:(n- 1)\100;
y=0:100:(m- 1)\100;
subplot(2,1,1);
mazz=max(max(z))\3;
dip=40;dir=340;
meshz(z);
view(dir,dip);
axis([0 n- 1 0 m- 1 0 mazz]);
nn = (n- 1)\(n- 1);
mm =(m- 1)\(m- 1);
xl = x(n)- x(1);
yl = y(m)- y(1);
a=(x(2:n)- x(1:n- 1))/xl;
b=(x(n)*x(1:n- 1)- x(1)5*(2:n))/xl;
c=(y(2:m)- y(1:m- 1))/yl;
d=(y(m)*y(1:m- 1)- y(1)*y(2:m))/yl;
cz=z(1,1)+z(m,n)- z(1,n)- z(m,1);
cm=x(1)*y(1)+x(n)*y(m)- x(n)*y(1)- x(1)*y(m);
bz1=z(1,1)- z(1,n); 
bz2=x(1)*y(1)- x(n)*y(1);
bm=x(1)- x(n);
dz1=z(1,1)- z(m,1);
dz2=x(1)*y(1)- x(1)*y(m);
dm=y(1)- y(m);
dn=ones(m- 1,n- 1);
dn=dn*af;
cc=(z(1:m- 1,1:n- 1)- z(1:m- 1,2:n)- z(2:m,1:n- 1)+z(2:m,2:n)- dn*cz)/cm;
bb=(z(1:m- 1,1:n- 1)- z(1:m- 1,2:n)- dn5bz1- cc*bz2)/bm;
dd=(z(1:m- 1,1:n- 1)- z(2:m,1:n- 1)- dn5dz1- cc*dz2)/dm;
kk=z(2:m,2:n)- bb*x(n)- dd*y(m)- dn*z(m,n)- cc*x(n)*y(m);
for j=1:m- 1
    for j0=1:m
        yv=c(j)*y(j0)+d(j);
        jj=(j- 1)*(m- 1)+j0;
        for i=1:n- 1
            for i0=1:n
                ii=(i- 1)*(n- 1)+i0;
                xv=a(i)*x(i0)+b(i);
                zt=bb(j,i)*x(i0)+dd(j,i)*y(j0);
                zz(jj,ii)=zt+cc(j,i)5x(i0)*y(j0)+dn(j,i)*z(j0,i0)+kk(j,i);
            end
        end
    end
end
mm=(m- 1)*(m- 1)+1;
nn=(n- 1)*(n- 1)+1;
subplot(2,1,2);
meshz(zz);
view(dir,dip);
axis([0 nn- 1 0 mm- 1 0 mazz]);
axis off
spp=['fwd=fopen(' '''' sp '''' ',' '''' 'w' '''' ')'];
eval ([spp]);
for j=1:mm
    for i=1:nn
        fprintf(fwd,'%8.4f',zz(j,i));
    end%
    fprintf(fwd,'\n');%
end%
fclose(fwd);%
end