我试图使用numpy产生2D柏林噪声，但我得到的不是平滑的东西：
我破碎的柏林噪音，到处都是丑陋的方块
当然，我在某个地方混淆了我的维度，可能是当我把四个梯度组合在一起的时候。。。但我找不到它，我的大脑正在融化。有人能帮我找出问题吗？
无论如何，这是代码：
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
def perlin(x,y,seed=0):
# permutation table
np.random.seed(seed)
p = np.arange(256,dtype=int)
np.random.shuffle(p)
p = np.stack([p,p]).flatten()
# coordinates of the first corner
xi = x.astype(int)
yi = y.astype(int)
# internal coordinates
xf = x - xi
yf = y - yi
# fade factors
u = fade(xf)
v = fade(yf)
# noise components
n00 = gradient(p[p[xi]+yi],xf,yf)
n01 = gradient(p[p[xi]+yi+1],xf,yf-1)
n11 = gradient(p[p[xi+1]+yi+1],xf-1,yf-1)
n10 = gradient(p[p[xi+1]+yi],xf-1,yf)
# combine noises
x1 = lerp(n00,n10,u)
x2 = lerp(n10,n11,u)
return lerp(x2,x1,v)
def lerp(a,b,x):
"linear interpolation"
return a + x * (b-a)
def fade(t):
"6t^5 - 15t^4 + 10t^3"
return 6 * t**5 - 15 * t**4 + 10 * t**3
def gradient(h,x,y):
"grad converts h to the right gradient vector and return the dot product with (x,y)"
vectors = np.array([[0,1],[0,-1],[1,0],[-1,0]])
g = vectors[h%4]
return g[:,:,0] * x + g[:,:,1] * y
lin = np.linspace(0,5,100,endpoint=False)
y,x = np.meshgrid(lin,lin)
plt.imshow(perlin(x,y,seed=0))
多亏了Paul Panzer和睡个好觉，它现在可以工作了。。。
import numpy as np
import matplotlib.pyplot as plt
def perlin(x, y, seed=0):
# permutation table
np.random.seed(seed)
p = np.arange(256, dtype=int)
np.random.shuffle(p)
p = np.stack([p, p]).flatten()
# coordinates of the top-left
xi, yi = x.astype(int), y.astype(int)
# internal coordinates
xf, yf = x - xi, y - yi
# fade factors
u, v = fade(xf), fade(yf)
# noise components
n00 = gradient(p[p[xi] + yi], xf, yf)
n01 = gradient(p[p[xi] + yi + 1], xf, yf - 1)
n11 = gradient(p[p[xi + 1] + yi + 1], xf - 1, yf - 1)
n10 = gradient(p[p[xi + 1] + yi], xf - 1, yf)
# combine noises
x1 = lerp(n00, n10, u)
x2 = lerp(n01, n11, u) # FIX1: I was using n10 instead of n01
return lerp(x1, x2, v) # FIX2: I also had to reverse x1 and x2 here
def lerp(a, b, x):
"linear interpolation"
return a + x * (b - a)
def fade(t):
"6t^5 - 15t^4 + 10t^3"
return 6 * t**5 - 15 * t**4 + 10 * t**3
def gradient(h, x, y):
"grad converts h to the right gradient vector and return the dot product with (x,y)"
vectors = np.array([[0, 1], [0, -1], [1, 0], [-1, 0]])
g = vectors[h % 4]
return g[:, :, 0] * x + g[:, :, 1] * y
# EDIT : generating noise at multiple frequencies and adding them up
p = np.zeros((100,100))
for i in range(4):
freq = 2**i
lin = np.linspace(0, freq, 100, endpoint=False)
x, y = np.meshgrid(lin, lin) # FIX3: I thought I had to invert x and y here but it was a mistake
p = perlin(x, y, seed=87) / freq + p
plt.imshow(p, origin='upper')
编辑(2023)：这篇文章似乎很受欢迎，所以我重新审视了一下。以前，代码在给定种子的情况下以一个频率生成噪声。
在这个新版本中，我添加了不同频率和振幅的噪音。这里，我使用的是频率[1,2,4,8]，振幅是频率的倒数。这样，低频率定义了整体形状，而高频率则增加了细节。
看到你最近在本月对文章进行了编辑，我想分享一个我重构的代码版本，以更好地理解分形噪声生成的基本原理：我修改了代码以删除任何显式的2D相关代码，这样它就能够为任何维度生成噪声。此外，我还添加了噪声的1D和3D切片的可视化，以及用于进一步优化代码的探查器。使用profiler，我添加了一些调整，以稍微加快噪音生成功能。然而，最大的改进可以是不多次重新计算同一网格位置的嵌套排列表查找，但到目前为止，我还没有找到解决这个问题的优雅方法。
为了回答@Tisklon的问题，此版本的允许您定义所需的任何梯度向量，只需更改get_gradients函数即可返回任何您喜欢的梯度向量。默认情况下，它返回所选维度中的所有对角向量。
#****************************************************************************************************
# Imports
#****************************************************************************************************
#--- Built-Ins ---
import itertools
from functools import cache
#--- Mathematics ---
import numpy as np
#****************************************************************************************************
# Noise Generation
#****************************************************************************************************
class NoiseGenerator():
#================================================================================
# Initialization
#================================================================================
def __init__(self,seed=0,nr_dimensions=2):
self.SEED = int(seed)
self.NR_DIMENSIONS = int(nr_dimensions)
self.fade_f = smootherstep
self.compute_constants()
def compute_constants(self):
self.PERMUTATION_TABLE = get_permutation_table(self.SEED)
self.CORNERS = get_corners (self.NR_DIMENSIONS)
self.GRADIENTS = get_gradients(self.NR_DIMENSIONS)
# Extend memory, to avoid '%' operation when retrieving gradient indices!
self.NR_GRADIENTS = self.GRADIENTS.shape[0]
GRADIENT_MULTIPLIER = int(np.ceil(self.PERMUTATION_TABLE.shape[0]/self.NR_GRADIENTS))
self.GRADIENTS_EXTENDED = np.vstack([self.GRADIENTS]*GRADIENT_MULTIPLIER)
#================================================================================
# Generation
#================================================================================
def fractal_noise(self,pos,octaves=8):
noise = np.zeros(pos.shape[:-1])
for i in range(octaves):
freq = 2**i
amp = 1/freq
noise+= self.perlin_noise(pos*freq) * amp
return noise
def perlin_noise(self,pos):
pos_i = pos.astype(int) # Grid coordinates
pos_f = pos - pos_i # Local fractional coordinates
gradients = {tuple(c):self.get_gradients (pos_i+c) for c in self.CORNERS} # Grid gradients # ToDo: Remove duplicate computation!
n = [self.dot(gradients[tuple(c)],pos_f-c) for c in self.CORNERS] # Noise components
pos_ff = self.fade_f(pos_f) # Fade positions
for i in range(self.NR_DIMENSIONS): # Interpolate noise
n = [lerp(n1,n2, pos_ff[self.filter_axis(i)]) for n1,n2 in zip(n[:len(n)//2],n[len(n)//2:])]
return n[0]
#================================================================================
# Support Functions
#================================================================================
def get_pos_grid(self,dim=512):
return np.moveaxis(np.mgrid[[slice(0,dim)]*self.NR_DIMENSIONS],0,self.NR_DIMENSIONS)/dim
def get_gradients(self,pos):
return self.GRADIENTS_EXTENDED[self.get_gradients_idx(pos)]
def get_gradients_idx(self,pos):
gradient_idx = pos[self.filter_axis(0)]
for i in range(1,self.NR_DIMENSIONS):
gradient_idx = self.PERMUTATION_TABLE[gradient_idx+pos[self.filter_axis(i)]]
return gradient_idx
def dot(self,a,b):
return np.sum([a[self.filter_axis(i)]*b[self.filter_axis(i)] for i in range(self.NR_DIMENSIONS)],axis=0)
def filter_axis(self,axis):
SLICE_ALL = [slice(None)]*self.NR_DIMENSIONS
return tuple(SLICE_ALL+[axis])
#================================================================================
# Support functions
#================================================================================
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Constants
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@cache
def get_permutation_table(seed=0,N=512):
np.random.seed(seed)
p = np.arange(N//2, dtype=int)
np.random.shuffle(p)
p = np.stack([p]*2).flatten()
return p
def _get_combinations(nr_dimensions,vs):
return np.array(list(itertools.product(*zip(*[[v]*nr_dimensions for v in vs]))))
@cache
def get_corners(nr_dimensions):
return _get_combinations(nr_dimensions,[0,1])
@cache
def get_gradients(nr_dimensions):
return _get_combinations(nr_dimensions,[-1,+1])
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Transitions
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def lerp(a, b, r):
return a + r * (b - a)
def smootherstep(t):
t2 = t**2 # Pre-compute square, for faster computation
return t2*t * (6*t2 - 15*t + 10)
#****************************************************************************************************
# Test Code
#****************************************************************************************************
if __name__=="__main__":
#--- Imports ---
import matplotlib.pyplot as plt
import time
#--- Settings ---
PROFILE = False
NR_DIMENSIONS = 3
DIM = 2**(9-NR_DIMENSIONS)
#--- Computation ---
t0 = time.time()
ng = NoiseGenerator(nr_dimensions=NR_DIMENSIONS)
pos = ng.get_pos_grid(DIM)
if PROFILE:
import cProfile
import pstats
cProfile.run("ng.fractal_noise(pos)","pstats.ps")
p = pstats.Stats("pstats.ps")
p.sort_stats(pstats.SortKey.TIME).print_stats(10)
quit()
else:
noise = ng.fractal_noise(pos)
dt = time.time()-t0
print(f"Noise generated in {dt:.1f} s")
#--- Visualization ---
def nd_slice(nd):
return tuple([slice(None)]*nd+[0]*(NR_DIMENSIONS-nd))
if NR_DIMENSIONS>=1:
plt.figure("1D")
plt.plot(np.arange(noise.shape[0]),noise[nd_slice(1)],color="gray")
if NR_DIMENSIONS>=2:
plt.figure("2D")
plt.imshow(noise[nd_slice(2)],cmap="gray")
if NR_DIMENSIONS>=3:
n = noise[nd_slice(3)]
n_norm = (n-np.min(n))/(np.max(n)-np.min(n))
colors = np.zeros(n.shape+(4,))
for i in range(4): colors[:,:,:,i] = n_norm
ax = plt.figure("3D").add_subplot(projection="3d")
ax.set(xlabel="x", ylabel="y", zlabel="z")
ax.voxels(
*np.indices(np.array(n.shape)+1), n,
facecolors=colors,
linewidth=0.5
)
plt.show()