MATLAB插值：柏林噪声插值

这是MATLAB插值拟合系列的第四期，这里附上前几期的链接：

第三期：MATLAB数值拟合：最小二乘多项式拟合

第二期：matlab插值：拉格朗日插值

第一期：[数值分析拟合]Matlab三次样条插值拟合数据

首先在这篇文章之前请允许我声明一下这篇文章确实是参考别人的算法得到的（算法最初的创造者是Ken Perlin）所以并不算我的原创内容（算是Ken Perlin的），先讲算法的实现过程，再上实际代码和演示。这一篇文章的话确实不能算原创，（虽然说里面的代码全部原创自己编写），（其实后半篇真的差不多是翻译[3]，只是尽可能让更多人听得懂）（像我这种人如果转载或翻译肯定要联系原作者的，我确实不知该标哪种类型，希望C站以后能多加几种文章类型）

（这篇文章就算是允许我仅仅代码部分标识原创而其他算法什么的都是别人的，讲解之类的都是附加的或翻译的内容,非作者本人原创）

先给参考文章哈：

一天多以来，为了弄柏林噪声我翻了B站和C站上面的许多东西，虽然也找到了一些资源但确实没找到什么能够把算法讲清楚的，虽然说维基上面有一些东西，但我也只是根据维基上面的实现了一维的插值，许多文章的话也往往只是有代码没方法

[1] 柏林噪声（Perlin Noise）但是这个主要原因我没有弄清楚里面的lerp之后它是怎么进行二维插值的

[2] 一篇文章搞懂柏林噪声算法，附代码讲解虽说这个确实相对有用，但到后来二维插值的时候它就讲的不是很多了（后面开始划水了）

不过，最终Lucky的我还是终于找到了文章[2]的参考文章也就是[3]，这篇文章讲的很详细，也是它最终让我弄懂了算法过程

[3] The Perlin noise math FAQ (作者是Matt Zucker)

下面进入正题——

柏林噪声是一种应用非常广泛的插值噪声算法，包括应用于生成云雾效果，火焰效果,熔岩等与地形（当然说到这里就不得不提到Minecraft地形生成其实也是可以用到这个算法的），包括应用于PhotoShop中的云彩渲染效果，也应用于Unity中的地形,波浪生成等等方面

一、一维柏林噪声插值的原理

我们先介绍两个平滑函数Smoothstep和Fade:

Smoothstep和Fade都是在[0,1]上进行平滑的函数，Smoothstep的表达式

$Smoothstep(x)= \begin{Bmatrix} 0& x\leqslant 0\\ 3x^2-2x^3 & 0<x<1 \\ 1& x\geqslant 1 \end{Bmatrix}$ 在[0,1]内图像如下：

而Fade函数则是一个更高阶的拟合，公式为

$Fade(x) = \begin{Bmatrix} 0& x\leqslant 0 \\ 6x^5-15x^4+10x^3 & 0<x<1\\ 1 & x\geqslant 1 \end{Bmatrix}$

（这些具体的可以参见维基里面的smoothstep，那里面也有更高阶的公式）

对于一维的柏林噪声插值，其原理是在每一段上使用Smoothstep或者Faded函数：

由于是在(0,1)区间内插值，那么我们一般来说将其限制在(0,1)范围内,所以我们一般取

$x_0=floor(x)$ ,其中，floor是向负无穷取整，这样就保证了要插的 $x-x_0$ 在（0,1）范围内,然后用公式：

$y = x_0+ Smoothstep(x-x_0)*(y_1-y_0)$

其中y1是第二个节点上面的函数值,(当然，这里Smoothstep函数也可以使用Fade函数来进行代替)

于是我改造了一维柏林噪声插值函数，让它能够插值任意的给定向量：

我把公式“篡改”为 $y = x_0+ Smoothstep(\frac{x-x_0}{x_1-x_0})*(y_1-y_0)$ ,这样就对于任意步长适用了。

二、一维柏林噪声代码实现

先给出我附加编写的函数（我一般把他们专门放在一个文件夹里并添加到路径下来调用）

1. len.m(计算向量长度（因为用length有点麻烦）)

(有的注释我用英文写的虽然写的不好可能还有语法错误吧)

%返回一维数组长度(行数或列数其中大的一个，必须是一位数组)
function [length]=len(A)s=size(A);if s(2)==1length=s(1); %取行endif s(1)==1length=s(2);endif s(1)~=1 && s(2)~=1disp('必须是单行向量或单列向量')returnendend

2. Smoothstep.m

function f = Smoothstep(x)
% define the function smoothstep that could be probably used in the Perlin noise
if x<=0f=0;
elseif x<1f=3*x^2-2*x^3;elsef=1;end
end
end

3.Fade.m

function f = Fade(x)
% define the function smoothstep that could be probably used in the Perlin noise
if x<=0f=0;
elseif x<1f=6*x^5-15*x^4+10*x^3;elsef=1;end
endend

4. (这个是主函数，其中如果输入两个参数默认使用Fade进行插值)，我的函数尽可能提高了兼容性，第三个参数既可以给插值函数@Fade/@Smoothstep也可以给插值步长dx（默认0.01）,最多四个参数都输，如果没有输出则绘制插值图像

function [x2,y2] = PerlinNoise_1(x,y,dx,Func)
% 柏林噪声一维插值函数
% 函数输入量只能为@Fade或者@Smoothstep
if nargin <=3Func = @Fade;   % use the function handle to give it a definitiondx = 0.01;
endif nargin ==3trytest = dx;   %test if the dx is validcatchdx = 0.01;endtrytest = Func(0.8); %test if the Func is validcatchFunc = @Fade;    end
endx2 = min(x):dx:max(x);
y2 = zeros(1,len(x2));j = 2;
% assure that the initial value should not be zero and(the initial value
% of two function are equal,so the first j should be 2)
% this denotes the x(j-1) which indicates the initial value of per step
for i= 1:len(x2)if x2(i)>x(j)j=j+1;endy2(i) = y(j-1) + Func((x2(i)-x(j-1))/(x(j)-x(j-1)))*(y(j)-y(j-1));  % attention ((x2(i)-x(j-1))/(x(j)-x(j-1))) is the ratio of every step
end
if nargout == 0plot(x2,y2,'-r','LineWidth',2)
end

然后使用测试代码：

x = 0:1.2:12;    % 这个是为了说明在步长不等于1时仍然适用的
y = cos(x);
PerlinNoise_1(x,y);

即可得到一维插值图像：

三、二维柏林噪声插值原理

下面讲二维柏林噪声插值，这些我都是参考自文章[3]的

对于一个二维的网格，我们仍然要进行插值，此时我们划分晶格之后，对晶格的每个点生成一个三维的随机向量（这个只要随机就好，不需要有过多要求），我们称之为梯度向量

当然，随机梯度向量更好的选择方法是从 $[[1,1,0];[-1,1,0];[1,-1,0];[-1,-1,0];[1,0,1];[-1,0,1];[1,0,-1]; [-1,0,-1];[0,1,1];[0,-1,1];[0,1,-1];[0,-1,-1]]$

这十二个向量（其实也对应着一个正方体中心指向其本身12条边中点的12个向量）中随机选取一个

所以我又编写了RandomVector来生成选取的随机向量:（就是随机生成一个下标并且选取）

function k = RandomVector()  % generate the vector we need
b_ = [[1,1,0];[-1,1,0];[1,-1,0];[-1,-1,0];[1,0,1];[-1,0,1];[1,0,-1];...[-1,0,-1];[0,1,1];[0,-1,1];[0,1,-1];[0,-1,-1]];
a_ = randRange(1,12.99999);
k = b_(fix(a_),:);
end

其中我又编写了randRange.m来生成指定范围内的随机数（其中要求用a,b来指定范围，Size用来指定生成指定范围内随机数矩阵的规模（不填默认返回一个随机数，比如填写[3,2]，那么返回的将是一个3*2的随机数矩阵，不过如果只填一个数比如Size填写4,则生成4*4的指定范围随机数方阵））

function X = randRange(a,b,Size)
% generate a matrix that cantains the random number in a specify range
% the third Agrument should be a vector
% if there is just one parament in the argument "size",it will create a
if nargin == 2Size = 1;
endX = ones(Size);
i=1;while i <= numel(X) % numel is used to get the number of the items of the matrix% for every item in the matrixX(i) = a+rand(1)*(b-a);i=i+1;
endend
% 当然其实这个我写的算法不是很好的，好的算法是直接用a + rand(Size)*(b-a),这样就不用for循环了

这样就可以生成随机的梯度向量了

我们在每个晶格点生成随机向量（包括边界）（容我用画图来说明）

于是放到三维空间里面，用matlab进行生成，结果如下：

对于晶格中的每一个插值点，定义四个距离向量d1,d2,d3,d4：

那么我们就有公式：

$\begin{Bmatrix} d1 = (x-x_0,y-y_0)\\ d2 = (x-x_1,y-y_0)\\ d3 = (x-x_0,y-y_1)\\ d4 = (x-x_1,y-y_1) \end{Bmatrix}$

于是分别定义每个点对插值点的影响为 $g\cdot d$ (其中，g，d分别为梯度向量和距离向量)（下面图来自文章[3]）

,来作为相应的影响值。（关于它为什么用这个做影响值我也没看懂,我觉的是以随机为原则选取的,s,t,u,v你可以把它们看做是在格点处的函数值）

那么，我们需要给s,t,u,v进行一系列加权，来确定最终中间插值点的z值，于是我们考虑先将s,t和u,v分别进行加权

$S_x = Smoothstep(x-x_0)$ （在代码中我使用的是Fade）

那么在x轴（s-t）方向上的分配来说，有

$a = s + S_x*(t-s)$ 这是针对于 s-t 一列的影响（加权后得到的）

$b = u + S_x(v-u)$ 这是 u-v 一列的影响（加权）

这样做相当于把两行分别进行插值

然后使用 $S_y = Smoothstep(y-y_0)$ , (y方向上的柏林噪声)

这样就能在y方向上对a ,b（两行）进行加权平均，权值为Sy

即： $z = a +S_y *(b-a)$

这样我们就能得到插值图像：

（照例先上代码再放图：(RandomVector我放在下面了,这个上面也单发过一次)

% Perlin noise interpolation in the second dimension
x = 0:1:4;
y = x;
[x1,y1] = meshgrid(x,y);
dx = 0.05;% ----------柏林噪声插值部分----------------------------------V = zeros(3,len(x),len(y));   % 生成随机向量矩阵
for i = 1:len(x)for j = 1:len(y)V(:,i,j) = RandomVector();end
endx2 = min(x):dx:max(x);
y2 = min(y):dx:max(y);%---------V(:,i,j)代表了在i,j坐标处的向量----------
Z2 = zeros(len(x2),len(y2));  %构造初始的插值矩阵for i = 1:len(x2)-1for j = 1:len(y2)-1  % 对中间的点进行插值% 使用floor向-inf取整g1 = V(:,floor(x2(i)+1),floor(y2(j)+1));  % 取第一个向量g2 = V(:,floor(x2(i)+2),floor(y2(j)+1));  % 取第二个向量g3 = V(:,floor(x2(i)+1),floor(y2(j)+2));  % 取第三个向量g4 = V(:,floor(x2(i)+2),floor(y2(j)+2));  % 注意:得到的向量是列向量%----------距离向量-------------------d1 = [x2(i)-floor(x2(i)),y2(j)-floor(y2(j)),0];d2 = [x2(i)-floor(x2(i)+1),y2(j)-floor(y2(j)),0];d3 = [x2(i)-floor(x2(i)),y2(j)-floor(y2(j)+1),0];d4 = [x2(i)-floor(x2(i)+1),y2(j)-floor(y2(j)+1),0];s = d1*g1 ;t = d2*g2 ;u = d3*g3 ;v = d4*g4 ;Sx = Fade(x2(i)-floor(x2(i)));
% in this case, the x2(i)-floor(x2(i))<1,so we can use this function to caculate ita = s + Sx*(t-s);b = u + Sx*(v-u);Sy = Fade(y2(j)-floor(y2(j)));Z2(i,j) = a + Sy*(b-a);%         hold on
%         plot3([floor(x2(i)),floor(x2(i))+g1(1)],[floor(y2(j)),floor(y2(j))+g1(2)],[0,g1(3)],'-r','LineWidth',2)end
end
subplot(1,2,1)
grid off
% colormap(hot)
surf(x2,y2,Z2)
% shading interp  %这两句是呈现熔岩效果的
subplot(1,2,2)
grid off
contour(x2,y2,Z2)%---------------------------------------------------function k = RandomVector()  % generate the vector we need
b_ = [[1,1,0];[-1,1,0];[1,-1,0];[-1,-1,0];[1,0,1];[-1,0,1];[1,0,-1];...[-1,0,-1];[0,1,1];[0,-1,1];[0,1,-1];[0,-1,-1]];
a_ = randRange(1,12.99999);
k = b_(fix(a_),:);   % fix向0取整
end

插值效果图:

也可以转一转

或者改改colormap来制造熔岩效果（这个代码我注释在代码里了）

不过这文章还没完哈

——为了说明s,t,u,v是随机影响(你可以把它看做某一点的随机函数值)，于是我把上面那一部分的s,t,u,v用3~4之间的随机数来代替，生成效果如下（避免数组溢出问题，我没有设置边值，所以4的地方函数值是0)

根据绘图结果显示函数值是在3~4之间的

MATLAB插值：柏林噪声插值相关推荐

MATLAB | MATLAB地形生成：矩形迭代法 · 傅里叶逆变换法 · 分形柏林噪声法
1:矩形迭代法这个非常简单,就是将矩阵的四个角分别定下初值,之后进行如下形式的迭代就好: [wxyz]⟶[ww+x2xw+y2w+x+y+z4x+z2yy+z2z]+noise.\begin{bma ...
matlab全域基函数,多项式函数插值：全域多项式插值（一）单项式基插值、拉格朗日插值、牛顿插值 [MATLAB]...
全域多项式插值指的是在整个插值区域内形成一个多项式函数作为插值函数.关于多项式插值的基本知识,见"计算基本理论". 在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值 ...
matlab 柏林噪声,游戏AI怎么写（一）——高级随机技术
写在前面最近在研究游戏AI,感觉要写一个不傻的AI并没有那么容易,于是开始研究AI到底该怎么写,有没有什么技巧和框架. 发现Web图书:Game AI Pro 刚开始读这一套书,我似乎意识到了把技巧 ...
利用MatLab对数据进行插值计算（分段插值和三次样条插值）
利用MatLab对数据进行插值计算分段线性插值三次样条插值例子分段线性插值应用的函数为: y=interp1(x0,y0,x)或y=interp1(x0,y0,x,'linear') 其中的 ...
解决龙格现象matlab,matlab实现Lagrange多项式插值观察龙格现象
Matlab进行Lagrange多项式插值拉格朗日插值法对函数y=1./(1+25*x.^2)在区间[-1,1]进行5次.10次.15次插值观察龙格现象主程序 1．拉格朗日 function [c ...
matlab插值法原理,【插值】插值方法原理详解
插值问题详解 1. 我在具体的应用(如数学建模竞赛)中,常常需要根据已知的函数点进行数据.模型的处理和分析,而通常情况下现有的数据是极少的,不足以支撑分析的进行,这时就需要使用一些数学的方法,&quo ...
清风数学建模学习笔记——应用matlab实现分段三次埃尔米特（Hermite）插值与三次样条插值
插值算法数模比赛中,常常需要根据已知的函数点进行数据.模型的处理和分析,而有时候现有的数据是极少的,不足以支撑分析的进行,这时就需要使用一些数学的方法,模拟产生一些新的但又比较靠谱的值来满足需求 ...
【matlab】三次埃尔米特插值与三次样条插值的实际应用代码
要求:完成下列这些数据的插值,并将结果保存到一个EXCEL表格中.要求至少选取两种插值方法,并对1号池中的这些指标做出插值后图像(显示在同一个图像中) Z.mat load Z.mat x=Z(1,: ...
Matlab实现线性插值、抛物插值、牛顿插值、拉格朗日插值、分段抛物插值、分段线性插值
目录线性插值原理流程图代码抛物插值原理流程图代码拉格朗日插值代码牛顿插值原理代码分段线性插值代码线性插值原理流程图单个点的线性插值代码 X=[0.2 0.4]; ...

MATLAB插值：柏林噪声插值

MATLAB插值：柏林噪声插值相关推荐

最新文章

热门文章