Deriving simplex noise

Abstract

The simplex noise algorithm, as described by Perlin in SIGGRAPH Course Notes (see references), has remained a bit mystical because the paper does not provide any derivations of the mathematical results. Furthermore, the provided source code is hard to understand because of its heavy optimizations and hard-coding of dimension-specific constants for $RR^3$ . To address some of these issues, Gustavson has written a paper to explain the algorithm more intuitively with source code that is easier to understand (see references). However, neither of the papers give a detailed derivation of the mathematical results in n-dimensional space. Such derivations are important in ensuring correctness and in establishing the claimed properties. In this paper we provide the following:

A derivation of the important results used in the simplex noise algorithm in n-dimensional space.
Observe that the support of the spherical attenuation function leaks into the neighboring simplices.
Give an attenuation function with a proper range.
Observe that with increasing dimensionality the algorithm generates increasing amounts of zero value. In the limit the space contains nothing else than the zero value.
Observe that the original algorithm is still exponential in complexity, although claimed otherwise.

Theory

Identifying simplices with vertex sets

Let $P = {p_0, p_1, ..., p_n} sub RR^n$ , where the points of $P$ are affinely independent. Then the convex hull of $P$ is a simplex. The convex hull of a $(k + 1)$ -sized subset of $P$ is a $k$ -simplex. In what follows, we identify a simplex with its vertices. A simplex is called regular if it holds that

$\exists d in RR: \forall i, j in [0, n]: i != j => |p_i - p_j| = d$

where the norm is Euclidean.

Simplicial partitioning of integer cubes

In the following, by a partitioning of $RR^n$ we mean a collection of sets which cover $RR^n$ and whose pair-wise intersection has Lebesgue measure zero. An integer cube is a set of the form $[h_1, h_1 + 1] xx ... xx [h_n, h_n + 1] sub RR^n$ , where $(h_1, ..., h_n) in ZZ^n$ . An integer cube is identified with its lower integer coordinates $(h_1, ..., h_n) in ZZ^n$ . The simplex noise algorithm assumes that $RR^n$ has been partitioned by (not necessarily regular) simplices. To find such a partitioning, first consider a partioning given by the set of by integer cubes. If $(h_1, ..., h_n) in ZZ^n$ is an integer cube, then clearly its vertices are given by ${h_1, h_1 + 1} xx ... xx {h_n, h_n + 1} sub ZZ^n$ . We would like to partition the cube with simplices such that each simplex contains as vertices $(h_1, ..., h_n)$ and $(h_1 + 1, ..., h_n + 1)$ . Those simplices are given by the set:

$S = {{h, h + e_{r_1}, ..., h + sum_{i = 1}^k e_{r_i}, ..., h + sum_{i = 1}^n e_{r_i}} in ZZ^{n + 1} : {r_1, ..., r_n} = \[1, n\] sub ZZ}$

where $e_i$ is the i:th standard basis vector. Because the number of permutations of $1, n$ is $n!$ , there are also $n!$ simplices in $S$ . This procedure gives a partitioning of $RR^n$ into simplices. However, these simplices are not regular nor close to being regular.

Scaling along a cube diagonal

Define the regularity of a simplex as the ratio of the length of its shortest edge and its longest edge. It is easily seen that the regularity of all the simplices in the simplicial partitioning of integer cubes is $1 / sqrt(n)$ . We would like to have a transformation class which would have a simple form, and then from that class find a transformation which maximizes regularity. Following Perlin, we shall consider a transformation which scales the space along the vector $p = [1, ..., 1]$ . Such a scaling is given by:

$f : RR^n -> RR^n: f(x) = x + (alpha - 1) (: x, p :) / (: p, p :) p = x + (alpha - 1) / n (sum_{i = 1}^n x_i) p$

where $alpha in RR$ , $alpha > 0$ , and $(: , :)$ is the dot product. Applying this scaling to the standard basis vectors $e_i$ gives:

$f(e_i) = e_i + (alpha - 1) / n p$

Edge lengths of the transformed simplex

A simplex from an integer cube partitioning is transformed into:

${f(h), f(h) + f(e_{r_1}), ..., f(h) + sum_{i = 1}^k f(e_{r_i}), ..., f(h) + sum_{i = 1}^n f(e_{r_i})}$

The vector between simplex vertices $k_1$ and $k_2$ , with $k_1 < k_2$ , is given by:

$D = (f(h) + sum_{i = 1}^{k_2} f(e_{r_i})) - (f(h) + sum_{j = 1}^{k_1} f(e_{r_j}))$

$= sum_{i = k_1 + 1}^{k_2} f(e_{r_i})$

$= sum_{i = k_1 + 1}^{k_2} (e_{r_i} + (alpha - 1) / n p)$

$= (k_2 - k_1) (alpha - 1) / n p + sum_{i = k_1 + 1}^{k_2} e_{r_i}$

Its squared length is given by:

$(: D :) = (: (k_2 - k_1) (alpha - 1) / n p + sum_{i = k_1 + 1}^{k_2} e_{r_i} :)$

$= (: (k_2 - k_1) (alpha - 1) / n p :) + (: sum_{i = k_1 + 1}^{k_2} e_{r_i} :) + 2 (: (k_2 - k_1) (alpha - 1) / n p, sum_{i = k_1 + 1}^{k_2} e_{r_i} :)$

$= (k_2 - k_1)^2 (alpha - 1)^2 / n^2 (: p :) + (k_2 - k_1) + 2 (k_2 - k_1) (alpha - 1) / n (: p, sum_{i = k_1 + 1}^{k_2} e_{r_i} :)$

$= (k_2 - k_1)^2 (alpha - 1)^2 / n + (k_2 - k_1) + 2 (k_2 - k_1) (alpha - 1) / n (k_2 - k_1)$

Denote $m = k_2 - k_1$ . Then:

$(: D :) = m^2 (alpha - 1)^2 / n + 2 m^2 (alpha - 1) / n + m$

$= m^2 ((alpha - 1)^2 + 2 (alpha - 1)) / n + m$

$= m^2 (alpha^2 - 1) / n + m$

Choosing $alpha$ to maximize regularity

We would like to choose $alpha$ such that regularity is maximized. Let

$g(m) = m^2 (alpha^2 - 1) / n + m$ .

Then

$g'(m) = 2 m (alpha^2 - 1) / n + 1$ .

The critical point of $g$ (assume $alpha != 1$ ) is found by equating its derivative to zero. Whether it is a maximum or a minimum depends on the value of $alpha$ .

$g'(m) = 0 <=> 2 m (alpha^2 - 1) / n + 1 = 0 <=> m = -n / (2(alpha^2 - 1)) <=> m = n / (2(1 - alpha^2))$

From here on the computation of regularity branches based on the position $m$ of the critical point. The following inequalities reveal the dependence of $m$ to $alpha$ :

$m >= 1 <=> n / (2(1 - alpha^2)) >= 1 <=> n / 2 >= 1 - alpha^2 <=> alpha^2 >= 1 - n / 2 <=> true$

$m <= n <=> n / (2(1 - alpha^2)) <= n <=> 1 / 2 <= 1 - alpha^2 <=> alpha^2 <= 1 / 2$

Thus, if it holds that $alpha^2 >= 1 / 2$ , then $g$ is increasing in $1, n$ and the squared regularity is given by:

$text(reg)_1^2 = g(1) / g(n) = ((alpha^2 - 1) / n + 1) / (n^2 (alpha^2 - 1) / n + n) = ((alpha^2 - 1) + n) / (n^2 (alpha^2 - 1) + n^2) = (alpha^2 - 1 + n) / (n^2 alpha^2)$

This is a decreasing function (w.r.t. $alpha$ ) when restricted to $alpha^2 >= 1 / 2$ , and thus it is maximized by $alpha^2 = 1 / 2$ :

$text(reg)_1^2 = ((-1 / 2) + n) / (n^2 / 2) = (2 n - 1) / n^2$

If it holds that $alpha^2 <= 1 / 2$ , then

$m_{text(min)} = { (1, text(, if ) alpha^2 >= 1 / (n + 1)), (n, text(, otherwise))]$

and

$m_{text(max)} = text(floor)(n / (2(1 - alpha^2)) + 1/2)$

If one picks $alpha^2 = 1 / (n + 1)$ , then

$text(reg)_2^2 = { (4 / (n + 2 + 1 / n), text(, if n is odd)), (4 / (n + 2), text(, if n is even))]$

For all $n >= 2$ it holds that $text(reg)_2 >= text(reg)_1$ . I suspect this is the value of $alpha$ which maximizes regularity. However, I haven’t managed to prove this yet. As a result of this scaling transformation all simplices have a minimum edge length of $d = sqrt(n / (n + 1))$ .

Deforming transformation

The transformation from the non-regular simplex partitioning to a more-regular simplex partitioning is linear and is given by the matrix $M in RR^{n xx n}$ :

$M = [f(e_1), ..., f(e_n)] = I + (1 / sqrt(n + 1) - 1) / n [p, ..., p] = I - (sqrt(n + 1) - 1) / (n sqrt(n + 1)) [p, ..., p]$

It can be shown that the inverse matrix of $M$ is given by:

$M^{-1} = I + (sqrt(n + 1) - 1) / n [p, ..., p]$

These are the transformations that seem to be used in the original algorithm. However, the problem with these transformations is that the minimum simplex edge length varies with dimension. This in turn implies that the feature size varies with dimension. Therefore, it is better to apply an additional uniform scaling, which expands distances but does not affect angles. Since the minimum length of the simplex edges is $d$ , the minimum edge length of 1 is obtained by multiplying with $1 / d$ :

$B = 1 / d M$

$B^{-1} = d M^{-1}$

Containing simplex

Assume the space is partitioned according to the more-regular simplicial partitioning as described above. Given a point $u in RR^n$ , we would like to obtain the vertices of that simplex which contains $u$ .

This is easily solved by transforming the vertices of the simplices to the integer grid. This works because containment is invariant to linear transforms. This transform was already given as the matrix $B^{-1}$ . It is worthwhile to expand the matrix-vector multiplication $B^{-1} u$ because of the simple structure of the matrix:

$x = [x_1, ..., x_n]^T = B^{-1} [u_1, ..., u_n]^T = d * ([u_1, ..., u_n]^T + (s * sum_{i = 1}^n u_i) * [1, ..., 1]^T)$

where $s = (sqrt(n + 1) - 1) / n$ . The simplicity of this transformation is not a coincidence: the simplex partitioning was specifically chosen (by Perlin) to yield this kind of transformation. The original algorithm does not have the $d$ factor.

Vertices of the containing simplex

Clearly the containing cube for a point $x$ is given by

$h = [h_1, ..., h_n]^T = text(floor)([x_1, ..., x_n]^T) = text(floor)(x)$ .

The simplices inside the cube can be identified with traversals in ${h_1, h_1 + 1} xx ... xx {h_n, h_n + 1} sub ZZ^n$ , such that the traversal begins from $[h_1, ..., h_n]$ and ends to $[h_1 + 1, ..., h_n + 1]$ , where one is only allowed to move along standard basis axes and no backing is allowed. It is easily seen that the number of different traversals is $n!$ , and that the path contains exactly $n + 1$ vertices.

Let $b = x - h$ , $b in \[0, 1\[^n sub RR^n$ . Further, order the components of $b$ into descending order ${b_{r_1}, ..., b_{r_n}}$ . Then the vertices of the containing simplex are given by ${h, h + e_{r_1}, ..., h + sum_{i = 1}^k e_{r_i}, ..., h + sum_{i = 1}^n e_{r_i}}$ . These (n + 1) vertices can be mapped back to the regular simplex space by using the transformation $B$ :

$u = [u_1, ..., u_n]^T = B [x_1, ..., x_n]^T = 1 / d * ([x_1, ..., x_n]^T - (q * sum_{i = 1}^n x_i) * [1, ..., 1]^T)$

where $q = (sqrt(n + 1) - 1) / (n sqrt(n + 1)) = s / sqrt(n + 1)$ .

Scalar fields

Consider the more-regular simplex space again. Each simplex vertex in this space is associated with a gradient vector. This is the gradient of the noise function at that point. The noise function itself is zero at each vertex. The gradient vector $g in RR^n$ of a given vertex $v in RR^n$ defines a scalar field $s$ centered on $v$ by:

$s : RR^n -> RR : s(x) = (: x, g :)$

This scalar field is attenuated by a spherically symmetric function $a$ centered on $v$ :

$a : RR^n -> RR : a(x) = {((1 - (: x :) / (r^2 d^2))^4, text(, if ) (: x :) <= r^2 d^2), (0, text(, otherwise))]$

where $r in \]0, 1\] sub RR$ is the percentage of the minimum simplex edge length $d$ after which $a$ must have been attenuated to zero. Because of the normalization of the simplex edge length, we have $d = 1$ and thus the expression for $a$ simplifies to:

$a : RR^n -> RR : a(x) = {((1 - (: x :) / r^2)^4, text(, if ) (: x :) <= r^2), (0, text(, otherwise))]$

The product of the scalar functions $s$ and $a$ , centered on $v$ , is given by:

$e: RR^n -> RR : e(x) = s(x) * a(x)$

Choosing $r$

The $r$ must be chosen appropriately: the support of the function $e$ must be such that it is contained in the union of those simplices which are connected to the given vertex. If this were not the case, then it would not be sufficient to evaluate only $(n + 1)$ surrounding vertices for a given point. Out of the possible values for $r$ we want the largest one. This is given by $r$ such that $r d = h$ , where $h$ is the height of a regular simplex of the minimum edge length $d$ :

$h = d sqrt((3 / 4)^{n - 1})$

Thus

$r = sqrt((3 / 4)^{n - 1})$

and

$r^2 = (3 / 4)^{n - 1}$

The idea is that the true height of the simplex must be greater than the height of a regular simplex with minimum edge length.

Simplex noise function

Let $V = {v_1, v_2, ...} sub RR^n$ be the set of all simplex vertices. The noise function is given by:

‘’text(noise)(x) : RR^n -> RR : text(noise)(x) = sum_{i = 1}^oo e(x - v_i)

Because $e$ has bounded support, this sum actually contains only $(n + 1)$ non-zero terms: these terms are those corresponding to the vertices of the simplex containing $x$ .

Problems with the original algorithm

Difference to the original attenuation function

The attenuation function we give differs from that given by Perlin for $RR^3$ :

$a_2 : RR^3 -> RR : a_2(x) = {(8(0.6 - (: x :))^4, text(, if ) (: x :) <= 0.6), (0, text(, otherwise))]$

However, this function is erroneous: its support leaks into neighboring simplices. This is because the squared height of the simplex in $RR^3$ is given by:

$h^2 = d^2 (3 / 4)^2 = (3 / 4)^3 ~= 0.42$

which is less than 0.6.

Diminishing support of $e$ in higher dimensions

From the given formulae it is evident that as dimension grows, more and more of the scalar field $e$ inside the simplex is zero. For example, in $RR^7$ , only 42% of the shortest edge of a simplex is covered by the scalar field of its either end-vertex. This means that 16% of the edge has $e$ zero, which implies that there is considerable amount of volume inside the simplex which has $e$ zero. In the limit, as dimension increases, $e$ is zero almost everywhere.

This is a fundamental restriction of the simplex noise algorithm: it isn’t suited for higher dimensions. However, an open question (to me) is: what is the highest dimension in which the results are still acceptable? Clearly the result in $RR^7$ and in higher dimensions is not acceptable. On the other hand, from practical tests it can be seen that the result in $RR^4$ (e.g. time-varying 3d-noise) is still good. The question then is if the maximum acceptable dimension is 4, 5, or 6?

Still of exponential complexity

One major feature of simplex noise was supposed to be that it does not have exponential dependence to dimension. However, the algorithm that Perlin uses to compute the pseudo-random gradient index seems (to me) to have exponential dependence to dimension. This is why we use a traditional table of permutation values instead (of exponential size), as given by the GradientField class.

Acknowledgements

‘Hagman’ and ‘Tonico’ from the newsgroup sci.math helped me with finding the transformed simplex edges with minimum and maximum length. The name of the thread was Minimize and maximize.

References

Real-Time Shading SIGGRAPH Course Notes, Chapter 2: Noise Hardware, Ken Perlin, 2001.

Simplex Noise Demystified, Stefan Gustavson, 2005.