基于matlab的分形插值程序（二维和三维都有）

% This program is for self-affine fractal interpolation surface clc; clear; % bds = input('输入文件名:', 's'); % af = input('输入纵向压缩比：'); input_file_name = 'n_3d'; % 文件名 portait_rate = 0.2 ; % 纵向压缩比 % 输入输出文件路径 input_path = ['./data/3d_data/' input_file_name]; output_path = ['./data/3d_data/' input_file_name '_output']; % 加载数据,并清除无用变量 eval(['load ' input_path]); eval(['z=' input_file_name ';']); eval(['clear ' input_file_name]); [row,columns] = size(z) % z的形状 % x,y的网格点坐标 x = [0:30/(columns-1):30] % 插值前x方向的网格点的坐标，有n个 y = [0:30/(row-1):30]; % 插值前y方向的网格点坐标，有m个 % x,y方向的网格数 x_number = (columns-1)*(columns-1); % 插值后x方向的网格数 y_number = (row-1)*(row-1); % 插值后y方向的网格数 x_den = x(columns)-x(1) % 关于x的分母 y_den = y(row)-y(1) % 关于y的分母 % 求x,y的公式的四个参数 a = (x(2:columns)-x(1:columns-1))/x_den; b = (x(columns)*x(1:columns-1)-x(1)*x(2:columns))/x_den; c = (y(2:row)-y(1:row-1))/y_den; d = (y(row)*y(1:row-1)-y(1)*y(2:row))/y_den; % 求插值后y坐标 for j=1:row-1 for j0=1:row jj = (j-1)*(row-1)+j0; yy(jj) = c(j)*y(j0)+d(j); end end % 求插值后x坐标 for i=1:columns-1 for i0=1:columns ii = (i-1)*(columns-1)+i0; xx(ii) = a(i)*x(i0)+b(i); end end % 求z的参数 cz = z(1,1) + z(row,columns) - z(1,columns) - z(row,1); cm = x(1) * y(1) + x(columns)*y(row) - x(columns)*y(1) - x(1)*y(row); bz1 = z(1,1) - z(1,columns); bz2 = x(1)*y(1) - x(columns)*y(1); bm = x(1) - x(columns); dz1 = z(1,1) - z(row,1); dz2 = x(1)*y(1) - x(1)*y(row); dm = y(1) - y(row); dn = ones(row-1,columns-1); dn = dn*portait_rate; % 公式参数 cc = (z(1:row-1,1:columns-1)-z(1:row-1,2:columns)-z(2:row,1:columns-1)+z(2:row,2:columns)-dn*cz)/cm; bb = (z(1:row-1,1:columns-1)-z(1:row-1,2:columns)-dn*bz1-cc*bz2)/bm; dd = (z(1:row-1,1:columns-1)-z(2:row,1:columns-1)-dn*dz1-cc*dz2)/dm; kk = z(2:row,2:columns)-bb*x(columns)-dd*y(row)-dn*z(row,columns)-cc*x(columns)*y(row); % 计算z for j=1:row-1 for j0=1:row jj=(j-1)*(row-1)+j0; for i=1:columns-1 for i0=1:columns ii = (i-1)*(columns-1)+i0; zt = bb(j,i)*x(i0)+dd(j,j)*y(j0); zz(jj,ii) = zt+cc(j,i)*x(i0)*y(j0)+dn(j,i)*z(j0,i0)+kk(j,i); end end end end % z1=500; % z2=1000; z1 = min(min(zz)) z2 = max(max(zz)) y_number=(row-1)*(row-1)+1; x_number=(columns-1)*(columns-1)+1; x1 = min(xx); x2 = max(xx); y1 = min(yy); y2 = max(yy); dip=40; dir=-20; [x,y]=meshgrid(xx,yy); surf(x,y,zz); view(dir,dip); axis([x1,x2,y1,y2,z1,z2]) view(dir,dip); grid off v=axis; hold on plot3([x1 x1],[y1 y1],[v(5) zz(1,1)],'k'); plot3([x2 x2],[y1 y1], [v(5) zz(1,x_number)],'k'); % 保存数据 spp = ['fwd=fopen(' '''' output_path '''' ',' '''' 'w' '''' ')']; eval([spp]); for j=1:y_number for i=1:x_number fprintf(fwd,'%8.2f',zz(j,i)); end fprintf(fwd,'\n'); end fclose(fwd);