# Hill Noise

A new approach to random noise generation with novel mathematical properties.

Estimated reading time: 13 mins

June 1, 2017Noise functions are a way to produce continuous pseudorandom values across one or more dimensions. The simplest example is a one-dimensional noise function that makes the y-value wobble up and down smoothly, but unpredictably, as you move along the x-axis. It’s also possible to produce 2-dimensional noise, where the amplitude of the noise depends on the x,y coordinates you’re sampling.

This can be used to create hilly terrain in a video game, or more generally, it’s used in computer graphics for procedurally generating textures or making bumpy surfaces. The original version of Minecraft used a 2D noise function to generate heights for terrain as a function of (x,y) coordinates, and later versions used a 3D noise function to generate terrain density as a function of (x,y,z) coordinates to produce overhangs and caves, as Notch describes here.

## Existing Noise Functions

There’s many different techniques for making noise. The simplest is called “Value Noise” and it works by generating a set of random points, and then interpolating between them. Most implementations store the random points in a fixed-width array, and points past the end of the array wrap around. Value noise looks like this:

```
class ValueNoise:
= 256
NUM_VALUES = 10
SPACING def __init__(self, random=random):
self.values = [random() for _ in range(ValueNoise.NUM_VALUES)]
def evaluate(self, x):
# Look up the nearest two values
= floor(x/ValueNoise.SPACING)
x1 = ceil(x/ValueNoise.SPACING)
x2 = self.values[x1 % ValueNoise.NUM_VALUES]
n1 = self.values[x2 % ValueNoise.NUM_VALUES]
n2 # Smoothly interpolate between them using cosine interpolation
= .5 - .5*cos((x/ValueNoise.SPACING % 1)*pi)
k return ((1-k)*n1 + k*n2)
```

More commonly, Perlin Noise and Simplex Noise are used. They are slightly more complicated algorithms that rely on similar principles: generate a large amount of random values, then look up and combine those values to get the value at a particular point. However, with Perlin and Simplex noise, the values are used as gradients (partial derivative vectors), rather than heights, which produces slightly better-looking noise.

### Fine Details

For many applications, it’s nice to have some extra detail. For example, if you’re using a noise function to make mountains, it’s good to have both “general trend” noise like “go upward for 50 miles, then go downward again”, but it’s also useful to have “local variation” noise like “there’s a bump on the road here”. Commonly, this is achieved by using multiple layers of noise. You start with high-amplitude, slow-changing noise (mountains) and then add smaller, quicker-changing noise (hills), and so on until you get to the smallest level of detail that you care about (bumps in the road). You can use the same noise function for each layer, if you scrunch it up by scaling the input and output of the function. The layers are commonly called “octaves”, and it looks like this:

```
class OctaveNoise:
def __init__(self, num_octaves, random=random):
self.noise = ValueNoise(random=random)
self.num_octaves = num_octaves
def evaluate(self, x):
= 0
noise for i in range(self.num_octaves):
= 0.5**i
size += size * self.noise.evaluate(x/size)
noise return noise/2
```

Unfortunately, most of the existing noise functions require pre-generating large amounts of random values that need to be interpolated between (commonly 256). If you keep sampling the noise function far enough in any one direction, the noise function will start repeating itself. Also, most noise functions don’t make any guarantees about the distribution of values. For example, Value Noise picks random points and interpolates between them, but most of the interpolated points are closer to 0.5 than to either 0 or 1, so the distribution of produced values is skewed towards the middle, and it’s quite rare to see the noise function get very high or very low, but quite common to see it pass through 0.5. Perlin Noise and Simplex Noise also suffer from these issues, though they do tend to produce better-looking results than Value Noise. For example, here’s what 1D Simplex Noise looks like:

Notice that the values produced by the Simplex Noise tend to be very heavily biased towards 0.5, and it’s very uncommon to see values close to 0 or 1. (The implementation here is based on the SimplexNoise1234 implementation used in the LÖVE game engine.) For many use cases, it might not matter that the distribution of values is not uniform. In my opinion though, it’s better to have a predictable distribution of values, which can be easily modified to get whatever distribution is desired for your use case.

## Hill Noise

As an alternative to address some of these issues, I present a new noise
function, which I call “Hill Noise”. It’s inspired by Fourier Decomposition: any periodic function can be broken
down into a sum of sine waves of different amplitudes, wavelengths, and offsets.
If you generalize the idea, you can create an arbitrarily complex and wobbly
function by summing together a bunch of sine waves of different amplitudes,
offsets, and wavelengths. If you sum together a bunch of sine waves, then the
pattern will only repeat itself after the least common multiple of the periods of the sine waves. For
example, the function
$f(x) = \sin(2\pi \frac{1}{2} x) + \sin(2\pi \frac{1}{3} x)$
has a period of
$6$.
But the function
$f(x) = \sin(2\pi \frac{100}{201} x) + \sin(2\pi \frac{100}{301} x)$
has a period of
$605.01$,
which is *much* larger because the two sine waves are slightly misaligned
and take a long time to resynchronize. For two sine waves with randomly chosen
periods, the resulting function typically has an incredibly large period. Even
if the periods aren’t randomly chosen, functions with a very long period can be
produced easily. For example:
$\sin(2\pi 0.91^0 x) + \sin(2\pi 0.91^1 x) + \sin(2\pi 0.91^2 x)$
has a period of 10,000.

```
class HillNoise:
def __init__(self, num_sines, random=random):
self.wavelengths = [random() for _ in range(num_sines)]
def evaluate(self, x):
= 0
noise for wavelength in self.wavelengths:
+= sin(x/wavelength)
noise return .5 + .5*noise/len(self.wavelengths)
```

The result of this is not terribly satisfying, because some of the sine waves are undesirably steep. If we’re trying to mimic the bumpiness of real life things, you typically don’t see bumps in the road that are as tall as a mountain, but as wide as a pebble. Generally, the shorter the length of a bump, the shorter the height of a bump. So, in order to emulate this, instead of $\sin(x/\texttt{wavelength})$, let’s use $\texttt{wavelength}\sin(x/\texttt{wavelength})$, which makes the amplitude directly proportional to the wavelength. This concept is similar to the Weierstrass Function, which has the sort of bumpiness we’re looking for.

```
class HillNoise:
def __init__(self, num_sines, random=random):
self.sizes = [random() for _ in range(num_sines)]
def evaluate(self, x):
= 0
noise for size in self.sizes:
+= size*sin(x/size)
noise return .5 + .5*noise/sum(self.sizes)
```

It’s a bit unfortunate that this noise function always evaluates to exactly 0.5 at x=0, so let’s introduce some random offsets. (This will be even more useful later.)

```
class HillNoise:
def __init__(self, num_sines, random=random):
self.sizes = [random() for _ in range(num_sines)]
self.offsets = [random()*2*pi for _ in range(num_sines)]
def evaluate(self, x):
= 0
noise for size,offset in zip(self.sizes, self.offsets):
+= size*sin(x/size + offset)
noise return .5 + .5*noise/sum(self.sizes)
```

Now, this function is already pretty useful, and it comes with built-in octaves with different levels of detail. But, ideally it would be nice if you could ask for different distributions of bumpiness. Maybe you don’t care about any noise whose amplitude is less than 0.1, or maybe you want a lot of small noise and a lot of big noise, but not much in between. So, let’s modify the function to take in an arbitrary set of wavelengths:

```
class HillNoise:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in range(len(sizes))]
def evaluate(self, x):
= 0
noise for size,offset in zip(self.sizes, self.offsets):
+= size*sin(x/size + offset)
noise return .5 + .5*noise/sum(self.sizes)
```

Because we added the random offsets earlier, we can now create two different noise functions with the same wavelengths, but which produce totally different noise! Nice!

So far, this is looking pretty good, but when we throw in a *whole
lot* of sine waves of similar amplitude, they kind of tend to average out
and we get a very bland noise function. According to the Central Limit Theorem, randomly sampling points on this
noise function will have an approximately normal distribution, and the standard
deviation will get smaller and smaller when more sine waves are added. In order
to counteract this, we need to use the cumulative distribution function of the
normal distribution. Unfortunately, it’s rather tricky to compute the cumulative
distribution function, but luckily we don’t need to be very exact, and Wikipedia
conveniently has an entire section on numerical approximations
of the normal CDF. Using a little bit of math black magic, it’s possible to
calculate the standard deviation of the sum of the sine waves, and plug it into
the CDF approximation.

```
class HillNoise:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in range(len(sizes))]
self.sigma = sqrt(sum((a/2)**2 for a in self.sizes))
def evaluate(self, x):
= 0
noise for size,offset in zip(self.sizes, self.offsets):
+= size*sin(x/size + offset)
noise # Approximate normal CDF:
/= sigma
noise return (0.5*(-1 if noise < 0 else 1)
*sqrt(1 - exp(-2/pi * noise*noise)) + 0.5)
```

The end result of this is noise that is still quite noisy, but now it’s (approximately) evenly distributed between 0 and 1:

With this evenly distributed noise function, it’s fairly easy to create noise with any precise distribution that you want using inverse transform sampling. For example, here’s a Weibull distribution (with adjustable parameters alpha and beta):

```
class WeibullHillNoise:
def __init__(self, sizes, alpha, beta, random=random):
self.hill_noise = HillNoise(sizes, random)
self.alpha, self.beta = alpha, beta
def evaluate(self, x):
= 1.0 - self.hill_noise.evaluate(x)
u return self.alpha * (-math.log(u)) ** (1.0/self.beta)
```

That’s quite powerful! And it’s not something that can be easily achieved with just any old noise function.

### 2D Hill Noise

Now, this is a lovely noise function, but it only works for 1 dimension of input. Let’s expand it to work in 2 dimensions! Instead of summing up sin(x) values, let’s sum up $\frac{\sin(x)+\sin(y)}{2}$ values.

```
class HillNoise2D:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in 2*range(len(sizes))]
self.sigma = sqrt(sum((a/2)**2 for a in self.sizes))
def evaluate(self, x, y):
= 0
noise for i,size in enumerate(self.sizes):
+= size/2*(sin(x/size + self.offsets[2*i])
noise + sin(y/size + self.offsets[2*i+1]))
# Approximate normal CDF:
/= 2*sigma
noise return (0.5*(-1 if noise < 0 else 1)
*sqrt(1 - exp(-2/pi * noise*noise)) + 0.5)
```

These will be axis-aligned, which is problematic because it will make the noise function much less random-looking, so for each 2D sine wave, let’s rotate its coordinates by a random angle.

```
class HillNoise2D:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in range(len(sizes))]
self.rotations = [random()*2*pi for _ in range(len(sizes))]
self.sigma = sqrt(sum((a/2)**2 for a in self.sizes))
def evaluate(self, x, y):
= 0
noise for i,size in enumerate(self.sizes):
# Rotate coordinates
= self.rotations[i]
rotation = x*cos(rotation) - y*sin(rotation)
u = -x*sin(rotation) - y*cos(rotation)
v += size/2*(sin(u/size + self.offsets[2*i])
noise + sin(v/size + self.offsets[2*i+1]))
# Approximate normal CDF:
/= 2*sigma
noise return (0.5*(-1 if noise < 0 else 1)
*sqrt(1 - exp(-2/pi * noise*noise)) + 0.5)
```

That looks pretty good, but sometimes the randomly generated rotations are nearly aligned for similar-sized sine waves, and it causes some unfortunate aliasing. Really, we don’t need the rotations to be random, but we want them to be fairly evenly distributed, and similarly sized sine waves should have not-parallel, not-perpendicular rotations. There’s a fairly clever trick for doing almost exactly this task using the Golden Ratio, that I learned about in this blog post about generating random colors programmatically. Let’s use that technique instead of purely random rotations.

```
= (sqrt(5)+1)/2
GOLDEN_RATIO class HillNoise2D:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in 2*range(len(sizes))]
self.sigma = sqrt(sum((a/2)**2 for a in self.sizes))
def evaluate(self, x, y):
= 0
noise for i,size in enumerate(self.sizes):
# Rotate coordinates
= (i*GOLDEN_RATIO % 1)*2*pi
rotation = x*cos(rotation) - y*sin(rotation)
u = -x*sin(rotation) - y*cos(rotation)
v += size/2*(sin(u/size + self.offsets[2*i])
noise + sin(v/size + self.offsets[2*i+1]))
# Approximate normal CDF:
/= 2*sigma
noise return (0.5*(-1 if noise < 0 else 1)
*sqrt(1 - exp(-2/pi * noise*noise)) + 0.5)
```

Here’s what the finished product looks like:

### 3D Hill Noise

Great! Now what about 3D noise? Well… it gets a fair bit more complicated. For one thing, it’s fairly hard to visualize 3D noise, because it has a different value at every point in 3D space. It’s a lot harder to notice if there are any artifacts. So instead, I’m going to visualize it as 2D noise and use time as the 3rd dimension to produce a noise function that wobbles randomly over time. The trick we used for evenly spacing the rotations last time isn’t going to work anymore, because now we need random polar coordinate rotations. With a lot of experimentation, and thanks to the excellent paper Golden Ratio Sequences For Low-Discrepancy Sampling, I managed to come up with a way to generate evenly spaced polar rotations that satisfy the same requirements we had for 2D noise. The end result is starting to get a bit hairy, but it looks like this:

```
= (sqrt(5)+1)/2
GOLDEN_RATIO class HillNoise3D:
def __init__(self, sizes, random=random):
self.sizes = sizes
self.offsets = [random()*2*pi for _ in range(len(sizes))]
self.sigma = sqrt(sum((a/2)**2 for a in self.sizes))
def evaluate(self, x, y):
= floor(log((resolution-1)*sqrt(5) + .5)/log(GOLDEN_RATIO))
fib_num = floor(.5 + (GOLDEN_RATIO**fib_num)/sqrt(5))
dec = floor(.5 + dec/GOLDEN_RATIO)
inc = 0
noise = 0
j for i in range(len(self.sizes)):
if j >= dec:
-= dec
j else:
+= inc
j if j >= len(self.sizes):
-= dec
j # Convert golden ratio sequence into polar coordinate unit vector
= ((i*GOLDEN_RATIO) % 1) * 2*pi
phi = acos(-1+2*((j*GOLDEN_RATIO) % 1))
theta # Make an orthonormal basis, where n1 is from polar phi/theta,
# n2 is roated 90 degrees along phi, and n3 is the
# cross product of the two
= sin(phi)*cos(theta), sin(phi)*sin(theta), cos(phi)
n1x,n1y,n1z = cos(phi)*cos(theta), cos(phi)*sin(theta), -sin(phi)
n2x,n2y,n2z # Cross product
= (n1y*n2z - n1z*n2y,
n3x,n3y,n3z *n2x - n1x*n2z,
n1z*n2y - n1y*n2x)
n1x# Convert pos from x/y/z coordinates to n1/n2/n3 coordinates
= n1x*x + n1y*y + n1z*z
u = n2x*x + n2y*y + n2z*z
v = n3x*x + n3y*y + n3z*z
w # Pull the amplitude from the shuffled array index ("j"), not "i",
# otherwise neighboring unit vectors will have similar amplitudes!
= self.sizes[j]
size # Noise is the average of cosine of distance along
# each axis, shifted by offsets and scaled by amplitude.
+= size/3*(cos(u/size + self.offsets[3*i])
noise + cos(v/size + self.offsets[3*i+1])
+ cos(w/size + self.offsets[3*i+2]))
# Approximate normal CDF:
/= 3*sigma
noise return (0.5*(-1 if noise < 0 else 1)
*sqrt(1 - exp(-2/pi * noise*noise)) + 0.5)
```

Whew! That’s a heck of a function! But look at the great results!

(The sizes of the sine waves here are generated by $(\frac{i+0.5}{N+1})^{-\log(smoothness)}$, where $\texttt{smoothness}$ is an adjustable parameter.)

With all of these changes, the resulting Hill Noise algorithm:

- Will (practically) never repeat itself.
- Has highly configurable characteristics that are easy to understand (sine wave counts and size distributions).
- Is memory efficient and easy to implement in a shader, even without access to a pseudorandom number generator.
- Produces noise that is (nearly) uniformly distributed across the [0,1] interval.
- Can be modified to calculate exact surface normals/gradients.