Perlin Noise

柏林噪声算法常用于程序生成随机内容，在游戏、电影等领域应用广泛。

本文的目标是使用 Python 生成如下的图像，由这个图像出发可以实现很多效果^[1]。

虽然柏林噪声算法有一个改进版本^[2]，但本文仍是基于旧原理的实现。

很多时候，我们并不是想要这样的随机：

我们需要更平滑的随机，就像本文第一张图片那样。

39-1.jpg（第一张图片）是一张灰度图，可以将其每个像素的灰度值映射为数值，用二维坐标索引，就得到一个三维的曲面（如果你愿意忽略它是离散的的话）。显然 39-1.jpg 就像是杂乱的水面，而 39-2.png 就是钉床了（虽然由于钉子非常密集躺上去应该还不至于当场去世）。

Perlin 的做法是：

在网格上生成图片，每个小格子称为“晶格”。首先，在格点上生成一个随机的梯度向量。对于晶格内的点，作由四个晶格格点到自身的向量，这些向量与对应的四个梯度向量作点积，得到格点上的四个数值，那么对于此晶格内的点，噪声值就是这四个数值的双线性插值^[3]。此外，你可能需要一个 $fade$ 函数使图像更平滑。

But wait … , why it even works ?

一个可能的解释：

对此我也不是很清楚。虽然这大致和梯度的几何意义有关（想象一下那个三维曲面），但是还是有一些细节不能理解。等以后搞明白了再说。

OK，原理介绍完了，缝一个简陋的代码：

import numpy as np
import imageio

def generate_gradient_field(m, n):
    """
    np.random.rand()是NumPy库中的一个函数
    用于生成一个具有指定维度的随机数数组
    此例中会生成一个形状为(m, n, 2)的三维数组
    """
    gradient_field = np.random.rand(m, n, 2)*2-1
    return gradient_field

def myfade(x)->float:
    return x*x*x*(x*(x*6-15)+10)  # 可能是自作多情的优化

def lerp(t, fa, fb):  # 线性插值
    return (1-t) * fa + t * fb

def single_point_noise(x, y, dx, dy, gradient_field):
    """
    x, y 为晶格编号
    dx, dy 相对坐标
    """
    x0 = int(x)  # 梯度向量编号
    y0 = int(y)
    x1 = x0 + 1
    y1 = y0 + 1  # 笛卡尔坐标系的坐标记法

    # 计算格点处的噪声值
    gradient00 = np.dot(gradient_field[y0, x0], [dx, dy])
    gradient01 = np.dot(gradient_field[y1, x0], [dx, dy - 1])
    gradient10 = np.dot(gradient_field[y0, x1], [dx - 1, dy])
    gradient11 = np.dot(gradient_field[y1, x1], [dx - 1, dy - 1])

    u = myfade(dx)
    v = myfade(dy)

    # 双线性插值
    x_interpolation = lerp(u, gradient00, gradient10)
    y_interpolation = lerp(u, gradient01, gradient11)
    noise_value = lerp(v, x_interpolation, y_interpolation)
    
    return noise_value

def fill_block(x, y):  # 需要离散化处理
    for i in range(1, l+1):
        for j in range(1, l+1):
            mymap[y*l+j][x*l+i] = single_point_noise(x, y, (i-0.5)/l, (j-0.5)/l, ggf)

x=eval(input("输入晶格横向个数："))
y=eval(input("输入晶格纵向个数："))
l=eval(input("输入每个晶格的像素数(l*l): "))
mymap = np.zeros([y*l+5, x*l+5], dtype = float, order = 'C')
# 边界冗余5像素防止数组越界

ggf = generate_gradient_field(x+2, y+2)

for i in range(x):
    for j in range(y):
        fill_block(i, j)

# 输出时去除冗余像素
imageio.imwrite("D:\mypycode\perlinnoise\\testpic1.jpg", mymap[:y*l, :x*l])

再放一张输出结果：

顺带一提，这种噪音图会引起博主的轻微不适，原因不明。

我们假设，图片一共被划分为 $x\cdot y$ 个晶格，横向 $x$ 个，纵向 $y$ 个。每个晶格是 $l\cdot l$ 像素的正方形。

向量场格点从 $0$ 开始编号，横：$0,1,2,\cdots ,x$ ；纵：$0,1,2,\cdots ,y$ .

晶格从 $0$ 开始编号，$(0,0),(0,1),\cdots ,(0,x-1);\cdots$

对于每个晶格，处理过程是类似的：

for i in range(x):
    for j in range(y):
        fill_block(i, j)

每个晶格里面，需要填充每个像素的噪声值。如何计算这个噪声值，上文已经提过。

其中，single_point_noise函数是在连续意义下的单位正方形中计算点的噪声值，而在生成图片时，却是离散的像素点。这里我直接简单粗暴地使用像素点的几何中心坐标，传值到single_point_noise函数中。另外，由于是单位正方形，所以值要除以 $l$ .

在此基础上，有一个简单的拓展：如果在时间维度上插值，就可以制作出三维的情形，或是动态变化的二维柏林噪声——这两者是一回事。

import numpy as np
import imageio

def generate_gradient_field(m, n, p):
    """
    np.random.rand()是NumPy库中的一个函数
    用于生成一个具有指定维度的随机数数组
    此例中会生成一个形状为(m, n, p, 3)的四维数组
    """
    gradient_field = np.random.rand(m, n, p, 3)*2-1
    return gradient_field

def myfade(x)->float:
    return x*x*x*(x*(x*6-15)+10)  # 可能是画蛇添足的优化

def lerp(t, fa, fb):  # 线性插值
    return (1-t) * fa + t * fb

def single_point_noise(x, y, t, dx, dy, dt, gradient_field):
    x0 = int(x)  # 梯度向量编号
    y0 = int(y)
    t0 = int(t)
    x1 = x0 + 1
    y1 = y0 + 1  # 笛卡尔坐标系的坐标记法
    t1 = t0 + 1

    # 计算格点处的噪声值
    gradient000 = np.dot(gradient_field[y0, x0, t0], [dx, dy, dt])
    gradient010 = np.dot(gradient_field[y1, x0, t0], [dx, dy-1, dt])
    gradient100 = np.dot(gradient_field[y0, x1, t0], [dx-1, dy, dt])
    gradient110 = np.dot(gradient_field[y1, x1, t0], [dx-1, dy -1, dt])
    gradient001 = np.dot(gradient_field[y0, x0, t1], [dx, dy, dt-1])
    gradient011 = np.dot(gradient_field[y1, x0, t1], [dx, dy-1, dt-1])
    gradient101 = np.dot(gradient_field[y0, x1, t1], [dx-1, dy, dt-1])
    gradient111 = np.dot(gradient_field[y1, x1, t1], [dx-1, dy -1, dt-1])

    u = myfade(dx)
    v = myfade(dy)
    w = myfade(dt)

    # 三线性插值
    x_interpolation0 = lerp(u, gradient000, gradient100)
    y_interpolation0 = lerp(u, gradient010, gradient110)
    noise_value0 = lerp(v, x_interpolation0, y_interpolation0)
    x_interpolation1 = lerp(u, gradient001, gradient101)
    y_interpolation1 = lerp(u, gradient011, gradient111)
    noise_value1 = lerp(v, x_interpolation1, y_interpolation1)
    noise_value = lerp(w, noise_value0, noise_value1)
    
    return noise_value

def fill_block(x, y, t):  # 需要离散化处理
    for i in range(1, l+1):
        for j in range(1, l+1):
            for k in range(1, l+1):
                mymap[y*l+j][x*l+i][t*l+k] = single_point_noise(x, y, t, (i-0.5)/l, (j-0.5)/l, (k-0.5)/l, ggf)

x=eval(input("输入晶格横向个数："))
y=eval(input("输入晶格纵向个数："))
t=eval(input("输入时间晶格个数："))
l=eval(input("输入每个晶格的像素数(l*l): "))
mymap = np.zeros([y*l+5, x*l+5, t*l+5], dtype = float, order = 'C')
# 边界冗余5像素防止数组越界

ggf = generate_gradient_field(x+2, y+2, t+2)

for i in range(x):
    for j in range(y):
        for k in range(t):
            fill_block(i, j, k)

output_folder = "D:\mypycode\perlinnoise\\timeaxis\\"  # 输出文件夹的路径

for t in range(1, t*l+1):
    # 提取灰度图像
    gray_image = mymap[:, :, t]
    gray_image = gray_image[:y*l, :x*l]  # 去除冗余像素

    output_path = f"{output_folder}gray_image_{t}.png"  # 输出文件的路径和文件名
    imageio.imwrite(output_path, gray_image)

将得到的图片转为 GIF（博客中的 PythonNote 有提），最终输出结果如下：

这个动图看起来动静似乎不明显，这是因为设置了t=1 ，过大的话计算量太大。将 t 再稍微放大点，gif 的时间间隔缩短点，效果就会好很多。

在拓展维度时，代码的编写很大程度上是出于对某种统一形式的信仰——这代码就该这么写，而不需要考察其细节。在编写代码的时候，我甚至没有查找三线性插值是如何进行的。

最后的 tip : 上面代码的输出文件名可能不太好，因为在合成 gif 时，顺序可能是按严格的字典序排的。我是最后手动修改了文件名，反正也不多，就没有改代码了，因为很讨厌处理字符串。

记录一下最近看的感觉不错的电影，因为实在没必要开个新文章写这些：壳中少女、末代皇帝。另外在看我推，目前看着还行。

---

^[1] The Absurd Usefulness of Noise in Game Development

^[2] Ken Perlin’s SIGGRAPH 2002 paper: Improving Noise （已备份）

^[3] 维基百科 双线性插值 （已备份）