digital_map.zip_matlab 地图插值_地图 _地形平滑_地形数据_模拟电子地图

%%电子地图产生算法 %第一步：随机地图的产生 %第二步：特征地形的模拟 %第三步：地形数据的插值 %第四步：地形数据的平滑 %------------------------------------------------------------------------% %第一步：随机地图的产生 %应用常量的定义 Txc = 1; % X方向的抗相关距离 Tyc = 1; % Y方向的抗相关距离 Sdterr = 1; % 二维随机过程的均方差 Lx = 119; % X方向的长度 Ly = 119; % Y方向的长度 Terrmn = 120;%期望的地形高度的均值 %Y方向地图边界自递归序列产生 ay = exp(-1/Tyc); variance_y = (Sdterr)^2 * (1 - (ay)^2); Zy = sqrt(variance_y/2) * randn(1,Ly); %X方向地图边界自递归序列产生 ax = exp(-1/Txc); variance_x = (Sdterr)^2 * (1 - (ax)^2); Zx = sqrt(variance_x/2) * randn(1,Lx); %将X和Y方向的两个序列组合成为边界 Zb = zeros(120); for i = 1:1:119 Zb(1,i+1) = Zy(i); end for j = 1:1:119 Zb(j+1,1) = Zx(j); end %通用二维一阶离散自递归产生地图内部信息 a1 = exp(-1/Txc); %a1,a2,a3为在递归过程中的三个参量 a2 = exp(-1/Tyc); a3 = -a1 * a2; variance_w = (Sdterr)^2 * (1 - a1^2) * (1 - a2^2); W = wgn(120,120,variance_w); for x = 2:1:120 for y = 2:1:120 Zb(x,y) = a1 * Zb(x-1,y) + a2 * Zb(x,y-1) + a3 * Zb(x-1,y-1) + W(x,y); end end %叠加希望的地形高度的均值Terrmn Zb = Zb + Terrmn; %将随机地图显示出来 Xb = 0:1:119; Yb = 0:1:119; Zb = reshape(Zb,120,120); [Xb Yb] = meshgrid(Xb,Yb); %%surf(Xb,Yb,Zb); %------------------------------------------------------------------------% %第二步：特征地形的模拟 Zi = [15 25 35 22 17];%在基准地形基础上第i个山峰的高度 x0 = [50.0 70.0 90.0 110.0 63.0]; y0 = [40.0 100.0 50.0 70.0 63.0]; xx =[5.0 12.0 8.0 7.0 6.0]; yy=[24.0 22.0 19.0 23.0 17.0]; syms x y z; Z0 = 120; for i = 1:5 Z0 = Z0 + Zi(i)*exp(-(x-x0(i))^2/(xx(i)^2)-(y-y0(i))^2/(yy(i)^2)); end X = 0:1:119; Y = 0:1:119; z_temp = subs(Z0,x,X); ZZ = subs(z_temp,y,Y); ZZ = reshape(ZZ,120,120); ZZ = ZZ + Zb; [XX YY] = meshgrid(X,Y); %%surf(XX,YY,ZZ); %------------------------------------------------------------------------% %第三步：地形数据的插值 xc = linspace(0,120,480);%将0——120之间插值为400个点 yc = linspace(0,120,480); [XC,YC] = meshgrid(xc,yc); ZC = interp2(X,Y,ZZ,XC,YC);%双线性插值算法函数 %surf(XC,YC,ZC); %------------------------------------------------------------------------% %第四步：地形数据的平滑 Lab = 0.25; angle_smax = pi/2.5; for xp = 1:1:479 for yp = 1:1:479 angle_ab = atan((ZC(xp,yp) - ZC(xp,yp+1))/Lab); if(angle_ab < -angle_smax) ZC(xp,yp) = ZC(xp,yp+1) - Lab * tan(angle_smax); end if(angle_ab > angle_smax) ZC(xp,yp+1) = ZC(xp,yp) - Lab * tan(angle_smax); end end end surf(XC,YC,ZC); %------------------------------------------------------------------------%