Gradient Noise

In the 1980s, Ken Perlin proposed an alternative noise algorithm that does not exhibit the artifacts of value noise. Perlin had been hired by the Mathematical Applications Group, Inc. (MAGI), to work on the movie Tron. The movie was a technical feat for the time, but Perlin said this of his experience:

Working on TRON was a blast, but on some level I was frustrated by the fact that everything looked machine-like. In fact, that became the "look" associated with MAGI in the wake of TRON. So I got to work trying to help us break out of the "machine-look ghetto."

Perlin published his noise algorithm in 1985. It was adopted by many artists in the motion picture industry, and Perlin went on to win an Academy Award for the algorithm in 1997.

Instead of generating random scalars at every lattice point, Perlin's algorithm generates random vectors. These vectors are considered to be tangent to an imaginary organic surface. Such vectors are called gradients, and Perlin's noise is therefore called gradient noise. The values themselves are not random; they are computed to fit the shape of the gradients.

TODO: image of 1D noise

Let's implement our own Perlin 2D noise generator. There are many tutorials on the web explaining Perlin noise. There are several reasons for their abundance: the algorithm is sophisticated enough to need explanation, Perlin himself published several different versions of the algorithm, and many explanations focus on optimizations. Our discussion here will stick to a very simple and unoptimized implementation of the algorithm.

Random Gradients

First we need a grid of random gradients pointing in all directions. We could generate a random gradient with code by choosing a random point within a unit square, like this:

function randomGradient
  x = random(-1, 1)
  y = random(-1, 1)
  return new Vector2(x, y).normalize

function randomGradient
  x = random(-1, 1)
  y = random(-1, 1)
  return new Vector2(x, y).normalize

There's a problem with this approach, however. We are more likely to choose a point on the square's diagonal than in a cardinal direction because there are more points on the diagonal. The gradients will not be uniformly distributed. To avoid this bias, we think of the gradient in polar coordinates—as a radius and angle—rather than in Cartesian coordinates. The radius is 1 since the gradient should have unit length. The angle is drawn randomly from $[0, 2\pi)$.

function randomGradient
  angle = random(0, 2 * pi)
  x = cos angle
  y = sin angle
  return new Vector2(x, y)

function randomGradient
  angle = random(0, 2 * pi)
  x = cos angle
  y = sin angle
  return new Vector2(x, y)

This pseudocode produces a uniformly distributed 2D array of gradients:

gradients = new grid of dimensions (width, height)
for each row r
  for each column c
    gradients[c, r] = randomGradient

gradients = new grid of dimensions (width, height)
for each row r
  for each column c
    gradients[c, r] = randomGradient

The result is a field of vectors that looks something like this:

We want to create this field of gradients just once and use it over and over again for all of our noise lookups. The exact resolution is not so important. It doesn't even need to match the resolution of the noise texture that we eventually generate. We'll make the gradients tile infinitely across the whole 2D plane.

Interpolation

The second ingredient to Perlin noise is a function that generates a random but coherent scalar value at an arbitrary position $p$. We need to find the fraction of $p$ within the cell, just as we did in the blerp routine in the Terrain class. Additionally, we need the coordinates of the four surrounding lattice points. All these values can be derived from the floor and ceiling of $p$:

function noise(p)
  floorP = floor of p
  ceilingP = floorP + [1, 1]
  fraction = p - floorP

function noise(p)
  floorP = floor of p
  ceilingP = floorP + [1, 1]
  fraction = p - floorP

This position may lie outside the gradient grid, so we wrap it to be in bounds.

function noise(p)
  // ...
  modFloorP = floorP % gradient dimensions
  modCeilingP = ceilingP % gradient dimensions

function noise(p)
  // ...
  modFloorP = floorP % gradient dimensions
  modCeilingP = ceilingP % gradient dimensions

The four neighboring gradients around $p$ shape the organic surface within the cell. We use the modded floor and ceiling grab the gradients at the four corners:

function noise(p)
  // ...
  swGradient = gradients[modFloorP.x, modFloorP.y]
  seGradient = gradients[modCeilingP.x, modFloorP.y]
  nwGradient = gradients[modFloorP.x, modCeilingP.y]
  neGradient = gradients[modCeilingP.x, modCeilingP.y]

function noise(p)
  // ...
  swGradient = gradients[modFloorP.x, modFloorP.y]
  seGradient = gradients[modCeilingP.x, modFloorP.y]
  nwGradient = gradients[modFloorP.x, modCeilingP.y]
  neGradient = gradients[modCeilingP.x, modCeilingP.y]

Somehow we need to turn $p$ and the four gradients into a random value that varies smoothly within the cell. Perlin's idea was to find vectors from the four corners to $p$:

function noise(p)
  // ...
  swVector = p - floorP
  seVector = p - [ceilingP.x, floorP.y]
  nwVector = p - [floorP.x, ceilingP.y]
  neVector = p - ceilingP

function noise(p)
  // ...
  swVector = p - floorP
  seVector = p - [ceilingP.x, floorP.y]
  nwVector = p - [floorP.x, ceilingP.y]
  neVector = p - ceilingP

And then dot each with the corresponding gradient vector:

function noise(p)
  // ...
  swDot = dot(swGradient, swVector)
  seDot = dot(seGradient, seVector)
  nwDot = dot(nwGradient, nwVector)
  neDot = dot(neGradient, neVector)

function noise(p)
  // ...
  swDot = dot(swGradient, swVector)
  seDot = dot(seGradient, seVector)
  nwDot = dot(nwGradient, nwVector)
  neDot = dot(neGradient, neVector)

These four dot products roughly describe how much $p$ aligns with the gradient vectors. We've got four numbers and want to reduce them down to one. Perlin suggested blerping them based on the fractional position of $p$:

function noise(p)
  // ...
  sDot = lerp(swDot, seDot, fraction.x)
  nDot = lerp(nwDot, neDot, fraction.x)
  dot = lerp(sDot, nDot, fraction.y)
  return dot

function noise(p)
  // ...
  sDot = lerp(swDot, seDot, fraction.x)
  nDot = lerp(nwDot, neDot, fraction.x)
  dot = lerp(sDot, nDot, fraction.y)
  return dot

The noise function produces only a single scalar value. To generate a complete field, we call it repeatedly across a grid:

function noiseField(width, height, scale)
  field = new field of dimensions (width, height)
  for each row r
    for each column c
      field[c, r] = noise((c, r) * scale)
  return field

function noiseField(width, height, scale)
  field = new field of dimensions (width, height)
  for each row r
    for each column c
      field[c, r] = noise((c, r) * scale)
  return field

The column and row indices are scaled, which zooms into the noise. The scale factor should be a number in $[0, 1]$.

With this much of the algorithm implemented, we are able to produce noise textures that look like this:

This isn't aesthetically pleasing. Let's fix it.

Normalize

The noise is currently too dark. That's because our noise values never reach 1. That would happen only if all the vectors were in alignment, and that's simply not possible. A workaround is to normalize the values to the range $[0, 1]$. We map the values from their native range $[\mathrm{min}, \mathrm{max}]$ to $[0, 1]$ with this arithmetic gauntlet:

$$\begin{align} [\mathrm{min}, \mathrm{max}] \xrightarrow{-~\mathrm{min}} [0, \mathrm{max} - \mathrm{min}] \xrightarrow{\div~(\mathrm{max} - \mathrm{min})} [0, 1] \end{align}$$

Enable normalization in this renderer to see the improvement:

Weight Smoothing

The normalized noise is more balanced, but the lightness reveals some strange discontinuities. These discontinuities appear between the cells of the gradient grid. Perlin identified their cause as the linear ramp of the fraction value of $p$. Linear ramps meet at a point, not a smooth transition. The fix is to run the fraction values through a smoothing function like this one:

$$\mathrm{smooth}(x) = 6x^5 - 15x^4 + 10x^3$$

This plot shows the raw fraction value in red and the smoothed fraction value in blue:

We apply the smoothing function to our fraction vector right before lerping:

function noise(p)
  // ...
  smoothFraction = [smooth(fraction.x), smooth(fraction.y)]
  sDot = lerp(swDot, seDot, smoothFraction.x)
  nDot = lerp(nwDot, neDot, smoothFraction.x)
  dot = lerp(sDot, nDot, smoothFraction.y)
  // ...

function noise(p)
  // ...
  smoothFraction = [smooth(fraction.x), smooth(fraction.y)]
  sDot = lerp(swDot, seDot, smoothFraction.x)
  nDot = lerp(nwDot, neDot, smoothFraction.x)
  dot = lerp(sDot, nDot, smoothFraction.y)
  // ...

Enable normalization and smoothing in this renderer to see a much nicer specimen of noise:

Octaves

With the normalization and smoothing enable, bump the resolution of the noise up to 200×200. The result is definitely better than white noise, but the texture lacks features. It looks worms crawled across some loose soil. We fix this by blending octaves of noise together. Just as before, big features must be assigned more weight and fine details less weight.

How do we determine these weights? Let's say the first octave has weight $\frac{1}{d}$, where $d$ is some denominator that we haven't figured out yet. Each successive weight is twice its predecessor, which gives us this geometric sequence:

$$\frac{1}{d}, \frac{2}{d}, \frac{4}{d}, \frac{8}{d} \ldots \frac{2^{c - 1}}{d}$$

We want the weights to sum to 1. This constraint means we can solve for the unknown denominator:

$$\begin{align} \frac{1}{d} + \frac{2}{d} + \frac{4}{d} + \frac{8}{d} + \ldots + \frac{2^{c - 1}}{d} &= 1 \\ \frac{1 + 2 + 4 + 8 + \ldots + 2^{c - 1}}{d} &= 1 \\ 1 + 2 + 4 + 8 + \ldots + 2^{c - 1} &= d \\ \end{align}$$

Denominator $d$ is a sum of a sequence of powers of 2. What do you notice about these sums of the first several sequences?

$$\begin{align} 1 &= 1 \\ 1 + 2 &= 3 \\ 1 + 2 + 4 &= 7 \\ 1 + 2 + 4 + 8 &= 15 \\ 1 + 2 + 4 + 8 + 16 &= 31 \\ 1 + 2 + 4 + 8 + 16 + 32 &= 63 \\ \end{align}$$

The sums are just one less than the next power of 2. In general, we compute $d$ with this shortcut:

$$\begin{align} d &= \sum_{i = 0}^{c - 1} 2^i \\ &= 2^c - 1 \\ \end{align}$$

Our fractalNoise function loops through the octaves, calls upon the regular noise function to produce each octave's random number, and computes a weighted sum of all the values:

fractalNoise(p, octaveCount) {
  sum = 0
  denominator = (1 << octaveCount) - 1
  for each octave i
    weight = (1 << i) / denominator
    sum += noise(p / (1 << i)) * weight
  return sum

fractalNoise(p, octaveCount) {
  sum = 0
  denominator = (1 << octaveCount) - 1
  for each octave i
    weight = (1 << i) / denominator
    sum += noise(p / (1 << i)) * weight
  return sum

The bit-shifting expression 1 << n is the same as the exponentiation 2 ^ n and is faster. Position $p$ is scaled down with each successive octave to ensure different positions are used. We update our noiseField function to call fractalNoise instead of noise.

Bump up the number of octaves in this final renderer to produce some very fine noise:

← Value Noise Heightmap Generation →

How to 3D