MATLAB | MATLAB地形生成：矩形迭代法 · 傅里叶逆变换法 · 分形柏林噪声法

1:矩形迭代法

这个非常简单，就是将矩阵的四个角分别定下初值，之后进行如下形式的迭代就好：
$$\left[\begin{array}{ccc}w& & x\\ & & \\ y& & z\end{array}\right]*\left[\begin{array}{ccc}w& \frac{w+x}{2}& x\\ \frac{w+y}{2}& \frac{w+x+y+z}{4}& \frac{x+z}{2}\\ y& \frac{y+z}{2}& z\end{array}\right]+\mathrm{noise}.$$
迭代过程示意图：

此方法的迭代法来自公众号Aki君，的如下推送(点击下方图片跳转链接)，此迭代法程序已经写的非常简练，并没有什么优化的必要，可以去自行查看：

function B = land(A)
% @author :aki
% @公众号 :Aki君

N = size(A,1);
d = (N+1)/2;   % dimension
level = log2(N-1);           % noise
scalef = 0.05*(2^(0.9*level));     

B = A;
B(d,d) = mean([A(1,1),A(1,N),A(N,1)]) + scalef*randn;   
B(1,d) = mean([A(1,1),A(1,N)]) + scalef*randn;           
B(d,1) = mean([A(1,1),A(N,1)]) + scalef*randn;            
B(d,N) = mean([A(1,N),A(N,N)]) + scalef*randn;
B(N,d) = mean([A(N,1),A(N,N)]) + scalef*randn;

if N>3
    B(1:d,1:d) = land(B(1:d,1:d));
    B(1:d,d:N) = land(B(1:d,d:N));
    B(d:N,1:d) = land(B(d:N,1:d));
    B(d:N,d:N) = land(B(d:N,d:N));
end
end

但是，迭代运算毕竟不是MATLAB的强项，当矩阵尺度逐渐增大时，迭代法运算速度会变慢的很快，毕竟创建进程会非常多，因此在这补充一个循环法，只有核心代码只有十行，不过可能确实不如迭代法容易理解：

function A=sqLand(A)
% @author :slandarer
% @公众号 :slandarer随笔

% 获取矩阵大小 N=2^k+1
N=size(A,1);
tcyc=log2(N-1);
% 循环数值填充
for i=tcyc:-1:1
    s=2^(i-1);scalef=0.05*(2^(0.9*i));
    A(s+1:2*s:N,1:s:N)=(A(1:2*s:N-1,1:s:N)+A(2*s+1:2*s:N,1:s:N))/2;
    A(1:s:N,s+1:2*s:N)=(A(1:s:N,1:2*s:N-1)+A(1:s:N,2*s+1:2*s:N))/2;
    A(s+1:2*s:N,1:s:N)=A(s+1:2*s:N,1:s:N)+scalef.*randn(2^(tcyc-i),2^(tcyc-i+1)+1);
    A(1:2*s:N,s+1:2*s:N)=A(1:2*s:N,s+1:2*s:N)+scalef.*randn(2^(tcyc-i)+1,2^(tcyc-i));
end
end

不过虽然不易理解但速度确实是快了很多，做个测试哈：

k=2^11+1;
A=zeros(k);
A([1 k],[1 k])=[1,1.25;1.1,2.0];

tic
B1=land(A);
toc

tic
B2=sqLand(A);
toc

历时 8.315744 秒。
历时 0.098721 秒。

调用函数绘图：

k=2^5+1;
A=zeros(k);
A([1 k],[1 k]) = [1,1.25;1.1,2.0];
% B=land(A);
B=sqLand(A);

colormap('copper')
strCell={
    'FontSize',12,'FontName','Cambria'};
Bisland=max(B,mean(B));
Bmax=max(max(Bisland));
Bmin=min(min(Bisland));

subplot(221),meshz(B)
axis([0 k 0 k min(min(B)) Bmax])
title('Default view',strCell{
    :})

subplot(222), meshz(Bisland)
axis([0 k 0 k Bmin Bmax])
title('Default view',strCell{
    :})

subplot(223), meshz(Bisland)
view([-75,40])
axis([0 k 0 k Bmin Bmax])
title('view([-75 40])',strCell{
    :})

subplot(224), meshz(Bisland)
view([240 55])
axis([0 k 0 k Bmin Bmax])
title('view([240 55])',strCell{
    :})

---

2:傅里叶逆变换法

来自公众号好玩的MATLAB这篇推送，微调了一下代码(点击下方图片跳转链接)，其实就是生成满足一定条件的随机系数矩阵，再做个二维傅里叶逆变换就好：

X=linspace(0,1,200)';
CL=(-cos(X*2*pi)+1).^.2;
r=(X-.5)'.^2+(X-.5).^2;
surf(X,X',abs(ifftn(exp(7i*rand(a))./r.^.9)).*(CL*CL')*30)

不过该算法比较适合生成单体景观，比较难做成多山峰的大型连续景观。

---

分形柏林噪声法

首先构造基本网格，再构造细分网格，基本网格上每一格点上都含有一个梯度向量：

细分网格上每一点都要与离其最近的四个基本网格上的梯度向量做内积。再将内积值乘以权重再相加即可获得噪声矩阵。

理论上这就是个二维插值，直接线性插值用上就完事：

但是，这样插值出来并不平滑，为了保证平滑性，我们需要将上述p,q值通过平滑的函数映射为其他值，常用的平滑函数有俩，一个是Smoothstep函数，此函数在x=0、x=0.5、x=1处分别等于0、0.5和1，0到1区间内整体也很平滑(f’(0)与f’(1)等于0)：

$$\mathrm{Smoothstep}(x)=\{\begin{array}{cc}0& x\le 0\\ 3x^2-2x^3& 0<x<1\\ 1& x\ge 0\end{array}$$

另一个则是Fade函数,不仅在x=0和x=1处的导数为0，而且导数的导数也为0，更加平滑:
$$\mathrm{Fade}(x)=\{\begin{array}{cc}0& x\le 0\\ 6x^5-15x^4+10x^3& 0<x<1\\ 1& x\ge 0\end{array}$$

对比一下：

基于上述描述，写出柏林噪声生成函数

function [Z,X,Y]=perlinNoise2(XLim,YLim,N,M,method,showNeg)
% @author :slandarer
% @公众号 :slandarer随笔

if nargin<6
    showNeg=false;
end

% 获取X,Y坐标点
[X,Y]=meshgrid(linspace(XLim(1),XLim(2),N+2),linspace(YLim(1),YLim(2),M+2));
X=X(2:end-1,2:end-1);Y=Y(2:end-1,2:end-1);

% 构造向量网格
if nargin<5,method=1;end
Msize=[diff(XLim)+1,diff(YLim)+1];
switch method
    case 1
        U=randi([1,2],Msize);
        V=(3-U)./sqrt(5)./2.*(randi([0,1],Msize)-0.5).*2;
        U=U./sqrt(5)./2.*(randi([0,1],Msize)-0.5).*2;
    case 2
        R=rand(Msize).*0.5;
        T=rand(Msize).*2.*pi;
        U=R.*cos(T);
        V=R.*sin(T);
end

Xl=floor(X);Xr=ceil(X);
Yd=floor(Y);Yu=ceil(Y);
% 获取网格左上角格点向量
Ilu=sub2ind(Msize,Xl-XLim(1)+1,Yu-YLim(1)+1);
Ulu=U(Ilu);Vlu=V(Ilu);
% 获取网格右上角格点向量
Iru=sub2ind(Msize,Xr-XLim(1)+1,Yu-YLim(1)+1);
Uru=U(Iru);Vru=V(Iru);
% 获取网格左下角格点向量
Ild=sub2ind(Msize,Xl-XLim(1)+1,Yd-YLim(1)+1);
Uld=U(Ild);Vld=V(Ild);
% 获取网格右下角格点向量
Ird=sub2ind(Msize,Xr-XLim(1)+1,Yd-YLim(1)+1);
Urd=U(Ird);Vrd=V(Ird);


fade=@(u)6.*u.^5-15.*u.^4+10.*u.^3;
% 做内积及映射
Zyu=(Ulu.*(X-Xl)+Vlu.*(Y-Yu)).*fade(Xr-X)+(Uru.*(X-Xr)+Vru.*(Y-Yu)).*fade(X-Xl);
Zyd=(Uld.*(X-Xl)+Vld.*(Y-Yd)).*fade(Xr-X)+(Urd.*(X-Xr)+Vrd.*(Y-Yd)).*fade(X-Xl);
Z=Zyu.*fade(Y-Yd)+Zyd.*fade(Yu-Y);

if showNeg
Z=Z+0.5;
end
Z(Z>1)=1;Z(Z<0)=0;
end

其中method取值为1或2，对应两种梯度向量生成方式，其一梯度向量将从(1,2),(2,1),(-1,2),(-2,1),(1,-2),(2,-1),(-1,-2),(-2,-1)这八个向量中随机选取，第二张方式则是圆内任意方向任意长度向量均可。

showNeg取值为true或false，意为是否将噪声曲面值增加0.5再舍去负值，以显示更多的凹陷部分。

调用函数代码：

% 显示单个柏林噪声
figure()
[Z,X,Y]=perlinNoise2([0 10],[0 10],100,100);
surf(X,Y,Z);

% 生成不同尺度柏林噪声
method=1;
[Z1,X1,Y1]=perlinNoise2([0 10],[0 10],120,120,method);
Z2=perlinNoise2([0 20],[0 20],120,120,method);
Z3=perlinNoise2([0 40],[0 40],120,120,method);
Z4=perlinNoise2([0 80],[0 80],120,120,method);
% 显示不同尺度柏林噪声
figure()
subplot(2,2,1)
pcolor(X1,Y1,Z1)
shading interp
subplot(2,2,2)
pcolor(X1,Y1,Z2)
shading interp
subplot(2,2,3)
pcolor(X1,Y1,Z3)
shading interp
subplot(2,2,4)
pcolor(X1,Y1,Z4)
shading interp

% 不同尺度柏林噪声叠加
figure()
p=0.3;% 不同尺度柏林噪声占比
surf(X1,Y1,Z1+Z2.*p+Z3.*p^2+Z4.*p^3);

% 模拟云雾
figure()
p=0.8;        % 不同尺度柏林噪声占比
method=1;     % 梯度向量生成方法
showNeg=true; % 是否保留负数部分
[Z1,X1,Y1]=perlinNoise2([0 10],[0 10],120,120,method,showNeg);
Z2=perlinNoise2([0 20],[0 20],120,120,method,showNeg);
Z3=perlinNoise2([0 40],[0 40],120,120,method,showNeg);
Z4=perlinNoise2([0 80],[0 80],120,120,method,showNeg);
pcolor(X1,Y1,Z1+Z2.*p+Z3.*p^2+Z4.*p^3)
shading interp
colormap("gray")

单个柏林噪声

显示不同尺度柏林噪声

不同尺度柏林噪声叠加

模拟云雾

---

完