1:矩形迭代法

这个非常简单，就是将矩阵的四个角分别定下初值，之后进行如下形式的迭代就好：
[wxyz]⟶[ww+x2xw+y2w+x+y+z4x+z2yy+z2z]+noise.\begin{bmatrix} w& &x\\ & & \\ y& &z \end{bmatrix}\longrightarrow \begin{bmatrix} 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{bmatrix}+\mathrm{noise}. ⎣⎡wyxz⎦⎤⟶⎣⎡w2w+yy2w+x4w+x+y+z2y+zx2x+zz⎦⎤+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>3B(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:1s=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);
toctic
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)：

Smoothstep(x)={0x≤03x2−2x30<x<11x≥0\mathrm{Smoothstep} (x)=\left\{\begin{array}{ccc} 0 &x\le 0\\ 3x^{2}-2 x^{3} & 0<\,x<1\\ 1 &x\ge 0 \end{array}\right. Smoothstep(x)=⎩⎨⎧03x2−2x31x≤00<x<1x≥0

另一个则是Fade函数,不仅在x=0和x=1处的导数为0，而且导数的导数也为0，更加平滑:
Fade(x)={0x≤06x5−15x4+10x30<x<11x≥0\mathrm{Fade} (x)=\left\{\begin{array}{ccc} 0 &x\le 0\\ 6 x^{5}-15 x^{4}+10 x^{3} & 0<\,x<1\\ 1 &x\ge 0 \end{array}\right. Fade(x)=⎩⎨⎧06x5−15x4+10x31x≤00<x<1x≥0

对比一下：

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

function [Z,X,Y]=perlinNoise2(XLim,YLim,N,M,method,showNeg)
% @author :slandarer
% @公众号 :slandarer随笔if nargin<6showNeg=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 methodcase 1U=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 2R=rand(Msize).*0.5;T=rand(Msize).*2.*pi;U=R.*cos(T);V=R.*sin(T);
endXl=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")

单个柏林噪声

显示不同尺度柏林噪声

不同尺度柏林噪声叠加

模拟云雾

完

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

【MATLAB】通信信号调制通用函数 — 傅里叶逆变换
目录傅里叶逆变换傅里叶逆变换 function [t,st] = F2T(f,sf) %This function calculate the time signal using ifft fun ...
matlab牛顿法求区间根程序,MATLAB用二分法、不动点迭代法及Newton迭代（切线）法求非线性方程的根...
一.实验原理二.实验步骤三.实验过程 1.(程序) (1)二分法:求在区间(1,2)之间的根,取 (a)bipart.m: function [x,m]=bipart(fun,a0,b0,to ...
重根的二阶迭代法matlab,MATLAB用二分法、不动点迭代法及Newton迭代（切线）法求非线性方程的根...
作者:凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/ 一.实验原理二.实验步骤三.实验过程 1.(程序) (1)二分法:求在区间(1,2)之间的根,取 ...
matlab中傅里叶反转亮度,Matlab傅里叶变换傅里叶逆变换-FFT-IFFT
<Matlab傅里叶变换傅里叶逆变换-FFT-IFFT>由会员分享,可在线阅读,更多相关<Matlab傅里叶变换傅里叶逆变换-FFT-IFFT(2页珍藏版)>请在人人文库网上搜 ...
【信号重构】经傅里叶逆变换(IFFT)后得到实数序列-含Matlab程序
我们处理的大部分数据都是实数序列.设有实数序列x(n),把x(n)经FFT转到频域中为X(k),然后在频域中进行处理,处理完后经IFFT变成实数序列.本实例将说明如何经IFFT后得到实数序列. 一.实 ...
matlab画梅花,基于Matlab图像素描生成算法究.doc
毕业文图像素描生成算法研究姓名院(系) 信息学院专业班级学号指导教师职称论文答辩日期年月日摘要分析比较图像处理提供参考.关键词: 目录 1 前言1 1.1 课题研究 ...
利用matlab实现复数域空间牛顿迭代法的分形图案展示（newton法）
利用matlab实现复数域空间牛顿迭代法的分形图案展示(newton法) 1 一维函数的牛顿迭代法 2 复平面的牛顿迭代法 2.1 简单方程结果 2.2 其它非线性方程结果本文首发于 matlab爱 ...
matlab迭代求解,[基于matlab平台的三种迭代法求解矩阵方程]matlab迭代法求方程的根...
数值分析第二次作业学院:电子工程学院基于matlab平台的三种迭代法求解矩阵方程组求解系数矩阵由16阶Hilbert方程组构成的线性方程组的解,其中右端项为[2877/851,3491/14 ...
matlab 生成连续信号,Matlab的连续信号生成及时频域分析
基于Matlab 的连续信号生成及时频域分析一.实验目的 1.通过实验使学生掌握matlab 表示信号的方法: 2.通过实验掌握基于matlab 的连续时间信号与系统的时频域分析方法. 二.实验要求 ...

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

1:矩形迭代法

2:傅里叶逆变换法

分形柏林噪声法

完

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

最新文章

热门文章