Value noise¶
A surface's variation, deterministic value noise in world space (solid noise), so it holds still as the camera moves and between steps.
At a point (x, y, z) in metres, the first octave's cells are scale metres wide: each corner (i, j, k) of the lattice of cells has a value in [0, 1), an integer hash of the corner and the octave (lusine.render.noise_corner), and the point's value is the cell's eight corners interpolated by smoothstep, u = f^2 (3 - 2 f) of the point's place f in its cell (along x, then z, then y), taken to [-1, 1]. Each next octave has cells half as wide and half the weight; the noise is the octaves' weighted sum over the weights' sum, octaves of them, at most 4 (0: none, 0). On the ground's plane y = 0 it is the plane's own value noise: a corner's hash with k = 0 is its hash of (i, j) alone, and a place of 0 across the cells' y takes the bottom corners' values as they are. Its gradient, lusine.render.value_noise_gradient, is the same sum's derivative in closed form: smoothstep's du/df = 6 f (1 - f), each octave's cells 2^-o scale wide, per metre. Source: Perlin, An Image Synthesizer (SIGGRAPH 1985), value noise by octaves and solid texture; Ebert et al., Texturing and Modeling: A Procedural Approach, 3rd ed., chapter 2 (lattice value noise). The hash is Wellons's "lowbias32" finaliser (hash-prospector, 2018) over the corner's coordinates times odd constants. Checked against a reference in numpy, which has its own tests.
- Source
- Perlin, An Image Synthesizer (SIGGRAPH 1985); Ebert et al., Texturing and Modeling: A Procedural Approach, 3rd ed., chapter 2; the hash's finaliser Wellons's lowbias32 (hash-prospector, 2018).
- Proof
- reference, invariant, limit.
value_noise.trait
namespace lusine.render;
// a whole number of cells: the cell holding x (floor), as an integer
integer noise_cell(float x)
{
return floor(x);
}
// a lattice corner's value in [0, 1): the hash of its cell coordinates and the octave, its top 24 bits over 2^24; a
// corner with cy = 0 hashes as the plane's corner (cx, cz)
float noise_corner(integer cx, integer cy, integer cz, integer octave)
{
let h0 = (bits(cx) * bits(2376512323)) ^ (bits(cz) * bits(3625334849)) ^ (bits(cy) * bits(3266489917))
^ (bits(octave) * bits(2654435769));
let h1 = h0 ^ (h0 >> 16);
let h2 = h1 * bits(2146121005);
let h3 = h2 ^ (h2 >> 15);
let h4 = h3 * bits(2221713035);
let h5 = h4 ^ (h4 >> 16);
return integer(h5 >> 8) * (1.0 / 16777216.0);
}
// one octave at x (in its cells), -1 to 1: the eight corners of x's cell smoothstep-interpolated, along x, then z,
// then y
float noise_octave(vector<float, 3> x, integer octave)
{
let ix = floor(x.x);
let iy = floor(x.y);
let iz = floor(x.z);
let fx = x.x - ix;
let fy = x.y - iy;
let fz = x.z - iz;
let cx = noise_cell(x = ix);
let cy = noise_cell(x = iy);
let cz = noise_cell(x = iz);
let ux = fx * fx * (3.0 - 2.0 * fx);
let uy = fy * fy * (3.0 - 2.0 * fy);
let uz = fz * fz * (3.0 - 2.0 * fz);
let a = noise_corner(cx = cx, cy = cy, cz = cz, octave = octave);
let b = noise_corner(cx = cx + 1, cy = cy, cz = cz, octave = octave);
let e = noise_corner(cx = cx, cy = cy, cz = cz + 1, octave = octave);
let g = noise_corner(cx = cx + 1, cy = cy, cz = cz + 1, octave = octave);
let a1 = noise_corner(cx = cx, cy = cy + 1, cz = cz, octave = octave);
let b1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz, octave = octave);
let e1 = noise_corner(cx = cx, cy = cy + 1, cz = cz + 1, octave = octave);
let g1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz + 1, octave = octave);
let lo = a + (b - a) * ux;
let hi = e + (g - e) * ux;
let lo1 = a1 + (b1 - a1) * ux;
let hi1 = e1 + (g1 - e1) * ux;
let v0 = lo + (hi - lo) * uz;
let v1 = lo1 + (hi1 - lo1) * uz;
return 2.0 * (v0 + (v1 - v0) * uy) - 1.0;
}
// one octave's gradient at x, per cell: the derivative of noise_octave along x, y and z
vector<float, 3> noise_octave_gradient(vector<float, 3> x, integer octave)
{
let ix = floor(x.x);
let iy = floor(x.y);
let iz = floor(x.z);
let fx = x.x - ix;
let fy = x.y - iy;
let fz = x.z - iz;
let cx = noise_cell(x = ix);
let cy = noise_cell(x = iy);
let cz = noise_cell(x = iz);
let ux = fx * fx * (3.0 - 2.0 * fx);
let uy = fy * fy * (3.0 - 2.0 * fy);
let uz = fz * fz * (3.0 - 2.0 * fz);
let dx = 6.0 * fx * (1.0 - fx);
let dy = 6.0 * fy * (1.0 - fy);
let dz = 6.0 * fz * (1.0 - fz);
let a = noise_corner(cx = cx, cy = cy, cz = cz, octave = octave);
let b = noise_corner(cx = cx + 1, cy = cy, cz = cz, octave = octave);
let e = noise_corner(cx = cx, cy = cy, cz = cz + 1, octave = octave);
let g = noise_corner(cx = cx + 1, cy = cy, cz = cz + 1, octave = octave);
let a1 = noise_corner(cx = cx, cy = cy + 1, cz = cz, octave = octave);
let b1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz, octave = octave);
let e1 = noise_corner(cx = cx, cy = cy + 1, cz = cz + 1, octave = octave);
let g1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz + 1, octave = octave);
let lo = a + (b - a) * ux;
let hi = e + (g - e) * ux;
let lo1 = a1 + (b1 - a1) * ux;
let hi1 = e1 + (g1 - e1) * ux;
let v0 = lo + (hi - lo) * uz;
let v1 = lo1 + (hi1 - lo1) * uz;
// along x: each edge's slope, then interpolated as the values are
let lo_x = (b - a) * dx;
let hi_x = (g - e) * dx;
let lo1_x = (b1 - a1) * dx;
let hi1_x = (g1 - e1) * dx;
let v0_x = lo_x + (hi_x - lo_x) * uz;
let v1_x = lo1_x + (hi1_x - lo1_x) * uz;
// along z and along y
let v0_z = (hi - lo) * dz;
let v1_z = (hi1 - lo1) * dz;
return 2.0 * vector<float, 3>(v0_x + (v1_x - v0_x) * uy, (v1 - v0) * dy, v0_z + (v1_z - v0_z) * uy);
}
float value_noise(vector<float<m>, 3> position, float<m> scale, integer octaves)
{
if (octaves <= 0) return 0.0;
let x = position / scale;
let n0 = noise_octave(x = x, octave = 0);
let n1 = octaves > 1 ? noise_octave(x = x * 2.0, octave = 1) : 0.0;
let n2 = octaves > 2 ? noise_octave(x = x * 4.0, octave = 2) : 0.0;
let n3 = octaves > 3 ? noise_octave(x = x * 8.0, octave = 3) : 0.0;
let total = ((n0 + 0.5 * n1) + 0.25 * n2) + 0.125 * n3;
let weights = octaves > 3 ? 1.875 : (octaves > 2 ? 1.75 : (octaves > 1 ? 1.5 : 1.0));
return total / weights;
}
// value_noise's gradient, per metre: octave o's weight 2^-o times its cells' 2^o per scale, so each octave's own
// gradient counts alike, over the weights' sum and the scale
vector<float<m^-1>, 3> value_noise_gradient(vector<float<m>, 3> position, float<m> scale, integer octaves)
{
let none_ = vector<float, 3>(0.0, 0.0, 0.0);
if (octaves <= 0) return none_ / scale;
let x = position / scale;
let g0 = noise_octave_gradient(x = x, octave = 0);
let g1 = octaves > 1 ? noise_octave_gradient(x = x * 2.0, octave = 1) : none_;
let g2 = octaves > 2 ? noise_octave_gradient(x = x * 4.0, octave = 2) : none_;
let g3 = octaves > 3 ? noise_octave_gradient(x = x * 8.0, octave = 3) : none_;
let weights = octaves > 3 ? 1.875 : (octaves > 2 ? 1.75 : (octaves > 1 ? 1.5 : 1.0));
return (((g0 + g1) + g2) + g3) / (weights * scale);
}
A surface's variation in world space, so it holds still as the camera moves: solid value noise, each corner of a lattice of cells a hash's value, interpolated by smoothstep, the octaves halving in cell and weight; on the ground's plane, y = 0, the plane's own noise of (x, z). Its gradient, value_noise_gradient, is the same sum's derivative in closed form, for a layer's bumps. One of the kernels lusine.render.table draws a table's shapes with, through a game's surface layers (lusine.render.surface).
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)
The kernel¶
The law as the library states it.
print(trait_run.source('render/value_noise.trait'))
// lusine.render.value_noise: a surface's variation, deterministic value noise in world space (solid noise), so it holds
// still as the camera moves and between steps. At a point (x, y, z) in metres, the first octave's cells are `scale`
// metres wide: each corner (i, j, k) of the lattice of cells has a value in [0, 1), an integer hash of the corner and
// the octave (lusine.render.noise_corner), and the point's value is the cell's eight corners interpolated by
// smoothstep, u = f^2 (3 - 2 f) of the point's place f in its cell (along x, then z, then y), taken to [-1, 1]. Each
// next octave has cells half as wide and half the weight; the noise is the octaves' weighted sum over the weights' sum,
// `octaves` of them, at most 4 (0: none, 0). On the ground's plane y = 0 it is the plane's own value noise: a corner's
// hash with k = 0 is its hash of (i, j) alone, and a place of 0 across the cells' y takes the bottom corners' values
// as they are.
// Its gradient, lusine.render.value_noise_gradient, is the same sum's derivative in closed form: smoothstep's du/df =
// 6 f (1 - f), each octave's cells 2^-o scale wide, per metre.
// Source: Perlin, An Image Synthesizer (SIGGRAPH 1985), value noise by octaves and solid texture; Ebert et al.,
// Texturing and Modeling: A Procedural Approach, 3rd ed., chapter 2 (lattice value noise). The hash is Wellons's
// "lowbias32" finaliser (hash-prospector, 2018) over the corner's coordinates times odd constants.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.render;
// a whole number of cells: the cell holding x (floor), as an integer
integer noise_cell(float x)
{
return floor(x);
}
// a lattice corner's value in [0, 1): the hash of its cell coordinates and the octave, its top 24 bits over 2^24; a
// corner with cy = 0 hashes as the plane's corner (cx, cz)
float noise_corner(integer cx, integer cy, integer cz, integer octave)
{
let h0 = (bits(cx) * bits) ^ (bits(cz) * bits) ^ (bits(cy) * bits)
^ (bits(octave) * bits);
let h1 = h0 ^ (h0 >> 16);
let h2 = h1 * bits;
let h3 = h2 ^ (h2 >> 15);
let h4 = h3 * bits;
let h5 = h4 ^ (h4 >> 16);
return integer(h5 >> 8) * (1.0 / 16777216.0);
}
// one octave at x (in its cells), -1 to 1: the eight corners of x's cell smoothstep-interpolated, along x, then z,
// then y
float noise_octave(vector<float, 3> x, integer octave)
{
let ix = floor(x.x);
let iy = floor(x.y);
let iz = floor(x.z);
let fx = x.x - ix;
let fy = x.y - iy;
let fz = x.z - iz;
let cx = noise_cell(x = ix);
let cy = noise_cell(x = iy);
let cz = noise_cell(x = iz);
let ux = fx * fx * (3.0 - 2.0 * fx);
let uy = fy * fy * (3.0 - 2.0 * fy);
let uz = fz * fz * (3.0 - 2.0 * fz);
let a = noise_corner(cx = cx, cy = cy, cz = cz, octave = octave);
let b = noise_corner(cx = cx + 1, cy = cy, cz = cz, octave = octave);
let e = noise_corner(cx = cx, cy = cy, cz = cz + 1, octave = octave);
let g = noise_corner(cx = cx + 1, cy = cy, cz = cz + 1, octave = octave);
let a1 = noise_corner(cx = cx, cy = cy + 1, cz = cz, octave = octave);
let b1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz, octave = octave);
let e1 = noise_corner(cx = cx, cy = cy + 1, cz = cz + 1, octave = octave);
let g1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz + 1, octave = octave);
let lo = a + (b - a) * ux;
let hi = e + (g - e) * ux;
let lo1 = a1 + (b1 - a1) * ux;
let hi1 = e1 + (g1 - e1) * ux;
let v0 = lo + (hi - lo) * uz;
let v1 = lo1 + (hi1 - lo1) * uz;
return 2.0 * (v0 + (v1 - v0) * uy) - 1.0;
}
// one octave's gradient at x, per cell: the derivative of noise_octave along x, y and z
vector<float, 3> noise_octave_gradient(vector<float, 3> x, integer octave)
{
let ix = floor(x.x);
let iy = floor(x.y);
let iz = floor(x.z);
let fx = x.x - ix;
let fy = x.y - iy;
let fz = x.z - iz;
let cx = noise_cell(x = ix);
let cy = noise_cell(x = iy);
let cz = noise_cell(x = iz);
let ux = fx * fx * (3.0 - 2.0 * fx);
let uy = fy * fy * (3.0 - 2.0 * fy);
let uz = fz * fz * (3.0 - 2.0 * fz);
let dx = 6.0 * fx * (1.0 - fx);
let dy = 6.0 * fy * (1.0 - fy);
let dz = 6.0 * fz * (1.0 - fz);
let a = noise_corner(cx = cx, cy = cy, cz = cz, octave = octave);
let b = noise_corner(cx = cx + 1, cy = cy, cz = cz, octave = octave);
let e = noise_corner(cx = cx, cy = cy, cz = cz + 1, octave = octave);
let g = noise_corner(cx = cx + 1, cy = cy, cz = cz + 1, octave = octave);
let a1 = noise_corner(cx = cx, cy = cy + 1, cz = cz, octave = octave);
let b1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz, octave = octave);
let e1 = noise_corner(cx = cx, cy = cy + 1, cz = cz + 1, octave = octave);
let g1 = noise_corner(cx = cx + 1, cy = cy + 1, cz = cz + 1, octave = octave);
let lo = a + (b - a) * ux;
let hi = e + (g - e) * ux;
let lo1 = a1 + (b1 - a1) * ux;
let hi1 = e1 + (g1 - e1) * ux;
let v0 = lo + (hi - lo) * uz;
let v1 = lo1 + (hi1 - lo1) * uz;
// along x: each edge's slope, then interpolated as the values are
let lo_x = (b - a) * dx;
let hi_x = (g - e) * dx;
let lo1_x = (b1 - a1) * dx;
let hi1_x = (g1 - e1) * dx;
let v0_x = lo_x + (hi_x - lo_x) * uz;
let v1_x = lo1_x + (hi1_x - lo1_x) * uz;
// along z and along y
let v0_z = (hi - lo) * dz;
let v1_z = (hi1 - lo1) * dz;
return 2.0 * vector<float, 3>(v0_x + (v1_x - v0_x) * uy, (v1 - v0) * dy, v0_z + (v1_z - v0_z) * uy);
}
float value_noise(vector<float<m>, 3> position, float<m> scale, integer octaves)
{
if (octaves <= 0) return 0.0;
let x = position / scale;
let n0 = noise_octave(x = x, octave = 0);
let n1 = octaves > 1 ? noise_octave(x = x * 2.0, octave = 1): 0.0;
let n2 = octaves > 2 ? noise_octave(x = x * 4.0, octave = 2): 0.0;
let n3 = octaves > 3 ? noise_octave(x = x * 8.0, octave = 3): 0.0;
let total = ((n0 + 0.5 * n1) + 0.25 * n2) + 0.125 * n3;
let weights = octaves > 3 ? 1.875: (octaves > 2 ? 1.75: (octaves > 1 ? 1.5: 1.0));
return total / weights;
}
// value_noise's gradient, per metre: octave o's weight 2^-o times its cells' 2^o per scale, so each octave's own
// gradient counts alike, over the weights' sum and the scale
vector<float<m^-1>, 3> value_noise_gradient(vector<float<m>, 3> position, float<m> scale, integer octaves)
{
let none_ = vector<float, 3>(0.0, 0.0, 0.0);
if (octaves <= 0) return none_ / scale;
let x = position / scale;
let g0 = noise_octave_gradient(x = x, octave = 0);
let g1 = octaves > 1 ? noise_octave_gradient(x = x * 2.0, octave = 1): none_;
let g2 = octaves > 2 ? noise_octave_gradient(x = x * 4.0, octave = 2): none_;
let g3 = octaves > 3 ? noise_octave_gradient(x = x * 8.0, octave = 3): none_;
let weights = octaves > 3 ? 1.875: (octaves > 2 ? 1.75: (octaves > 1 ? 1.5: 1.0));
return (((g0 + g1) + g2) + g3) / (weights * scale);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('render/value_noise.trait', 'value_noise'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float<m>, 3> position;
host float<m> scale;
host integer octaves;
float result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.render.value_noise(items.position, items.scale, items.octaves); }
from lusine.kernels.render.value_noise import value_noise as reference, value_noise_gradient as reference_gradient
Reference: a patch of ground¶
A 2 m patch of the plane at three octaves of 0.8 m: the kernel's noise within binary32's roundings of the float64 reference's, the hash the same integers in both.
def noise(position, scale=0.8, octaves=3):
return trait_run.law('render/value_noise.trait', 'value_noise', position=position, scale=scale, octaves=octaves)
def slope(position, scale=0.8, octaves=3):
return trait_run.law('render/value_noise.trait', 'value_noise_gradient', position=position, scale=scale,
octaves=octaves)
g = np.linspace(-1.0, 1.0, 51)
X, Z = np.meshgrid(g, g, indexing='ij')
pts = F32(np.stack([X.ravel(), np.zeros(X.size), Z.ravel()], axis=1)).astype(float)
got = noise(pts)
want = reference(pts, float(F32(0.8)), 3)
err = np.max(np.abs(got - want))
assert err <= 1e-5, err
fig, ax = plt.subplots(figsize=(4.5, 4))
ax.imshow(got.reshape(51, 51).T, origin='lower', extent=(-1, 1, -1, 1), cmap='gray')
ax.set_xlabel('x (m)'); ax.set_ylabel('z (m)')
plt.show()
print(f'{len(pts)} points within {err:.2g} of the reference')
2601 points within 1.9e-07 of the reference
Reference: solid, and its gradient¶
Points through a block of space, off the plane: the noise within binary32's roundings of the reference's, and its gradient (per metre) within the roundings its slopes carry, each octave's corners' differences over cells of a few centimetres.
block = F32(rng.uniform(-3.0, 3.0, (2000, 3))).astype(float)
err = np.max(np.abs(noise(block) - reference(block, float(F32(0.8)), 3)))
gerr = np.max(np.abs(slope(block) - reference_gradient(block, float(F32(0.8)), 3)))
scale = np.max(np.abs(reference_gradient(block, float(F32(0.8)), 3)))
assert err <= 1e-5 and gerr <= 1e-5 * scale, (err, gerr, scale)
print(f'noise within {err:.2g}; gradient within {gerr:.2g} m^-1 of slopes up to {scale:.2f} m^-1')
noise within 3.5e-07; gradient within 2.4e-06 m^-1 of slopes up to 3.38 m^-1
Invariant: bounded, centred and continuous¶
The noise stays in [-1, 1], averages about 0, and moves next to nothing for a micrometre's step (smoothstep's corners join with no jump); a step along the gradient moves it by the gradient's length times the step.
far = F32(rng.uniform(-50.0, 50.0, (2000, 3))).astype(float)
n0 = noise(far)
n1 = noise(far + 1e-4)
assert np.all(np.abs(n0) <= 1.0) and abs(float(np.mean(n0))) < 0.05
assert np.max(np.abs(n1 - n0)) < 2e-3
h = 1e-3
d = slope(far[:300])
fd = np.stack([(noise(far[:300] + h * e) - noise(far[:300] - h * e)) / (2 * h) for e in np.eye(3)], axis=1)
assert np.max(np.abs(d - fd)) < 0.05 * max(1.0, np.max(np.abs(d)))
print(f'mean {np.mean(n0):+.3f}, spread {np.std(n0):.3f}; the gradient the noise\'s slope')
mean +0.010, spread 0.237; the gradient the noise's slope
Limit: the plane's noise, the lattice's corners, and no octaves¶
On the plane y = 0 the solid noise is the plane's own value noise of (x, z), the table view's, exactly; at a corner of the first octave's lattice one octave is the corner's hash taken to [-1, 1]; no octaves is no noise, and no slope.
from lusine.kernels.render.value_noise import corner, _plane
flat = F32(rng.uniform(-20.0, 20.0, (2000, 3)) * np.array([1.0, 0.0, 1.0])).astype(float)
assert np.max(np.abs(noise(flat) - _plane(flat, float(F32(0.8)), 3))) <= 1e-5
corners = np.array([[0.0, 0.0, 0.0], [2.4, 0.8, -1.6], [-5.6, -1.6, 4.0]])
c = np.round(corners / 0.8).astype(int)
want = 2.0 * corner(c[:, 0], c[:, 1], c[:, 2], 0) - 1.0
assert np.allclose(noise(corners, 0.8, 1), want, rtol=0, atol=4 * U)
assert np.all(noise(corners, 0.8, 0) == 0.0) and np.all(slope(corners, 0.8, 0) == 0.0)
print('the plane\'s noise; a corner is its hash; no octaves, nothing')
the plane's noise; a corner is its hash; no octaves, nothing