n-dimensional interpolation

In the course of writing a perlin noise class which takes the dimension as a template parameter, I was in need of a function which interpolates in n dimensions. Generally speaking, interpolation is the process of “guessing” of a function value between two “known” function values, whereas extrapolation is guessing a new function value outside of the range of “known” values (see figure 4).

Figure 4: Interpolation and extrapolation

Figure 4: Interpolation and extrapolation

We first consider linear interpolation in one dimension. In this case, you have an \mathbb{R} vector space V of values – meaning the elements have the following two operators:

\alpha \cdot b: (\mathbb{R},V) \to V
a + b: (V,V) \to V

Those are the elements you want to interpolate. They could be colors, numbers, vectors, time stamps (!?) or whatever you can imagine. From this vector space, take two elements a_0, a_1 \in V and a scalar t \in \mathbb{R} and calculate

x = \textrm{interpolate}_1(t,a_0,a_1) = t \cdot a_0 + (1-t) \cdot a_1

Fig. 1: One-dimensional interpolation

Fig. 1: One-dimensional interpolation

In figure 1, this is depicted for V = \mathbb{R}^2. Now suppose you want to interpolate not on a line but on a square. In two dimensions, you have four points a_0,\ldots,a_3 instead of two (the number always doubles) and two scalars t_0,t_1 \in [0,1] (or one point t \in [0,1]^2), see figure 2 for an example.

Fig. 2: Two-dimensional interpolation

Fig. 2: Two-dimensional interpolation

To interpolate the point x between the four points, you have to do 2^n-1 = 3 interpolations, two for each line and then “inbetween” the lines.

x_0 = \textrm{interpolate}_1(t_0,a_0,a_1)
x_1 = \textrm{interpolate}_1(t_0,a_2,a_3)
x = \textrm{interpolate}_1(t_1,x_0,x_1)

And just for the fun of it, it’s also possible in three dimensions (so interpolation in a cube), where you have to interpolate first on the “front” side, then the “back” side and then inside the cube. See figure 3.

Figure 3: Three-dimensional interpolation

Figure 3: Three-dimensional interpolation

The corresponding interpolations are:

x_0 = \textrm{interpolate}_1(t_0,a_0,a_1)
x_1 = \textrm{interpolate}_1(t_0,a_2,a_3)
x_2 = \textrm{interpolate}_1(t_0,a_4,a_5)
x_3 = \textrm{interpolate}_1(t_0,a_6,a_7)
x_4 = \textrm{interpolate}_1(t_1,x_0,x_1)
x_5 = \textrm{interpolate}_1(t_1,x_2,x_3)
x = \textrm{interpolate}_1(t_2,x_4,x_5)

A general formula

Now as you can imagine, this can easily be generalized to n dimensions using some recursion formula. It took some time but I finally came up with one. Given a vector a = (a_0,\ldots,a_{2^n}) \in V^{2^n} of points and one “position” vector t \in [0,1]^n, define the interpolation in n dimensions as:

\textrm{interpolate}_n(a,t,i,f) := f(t_{n-1},\textrm{int}_0,\textrm{int}_1)
\textrm{interpolate}_1(a,t,i,f) := f(t_0,a_i,a_{i+1})


\textrm{int}_0 := \textrm{interpolate}_{n-1}(a,t,i,f)
\textrm{int}_1 := \textrm{interpolate}_{n-1}(a,t,i + 2^{n-1},f)

There are two new “unknowns” here, i and f. The first variable specifies the index to the vector a. It has to be 0 for the first function call. f specifies which “base” interpolation method to use. In the descriptions above we always assumed linear interpolation (so f = \textrm{interpolate}_1), but other methods are possible, like polynomial interpolation:

s(t) := 3t^2 - 2t^3
f(t,a_0,a_1) := s(t) \cdot a_0 + (1-s(t)) \cdot a_1

This gives a more smooth transition at the boundaries (in case you’re interpolating on a grid, see below), since the polynomial is smooth at 0 and 1, respectively. See figure 5.

Figure 5: Linear vs. polynomial interpolation

Figure 5: Linear vs. polynomial interpolation

To show that the formula makes sense, let’s apply it to the 2D case: We’ve got 4 values a_0,a_1,a_2,a_3 \in V, a point t=(t_0,t_1) and an interpolation function f. The first function call has to be with i=0. We get:

\textrm{interpolate}_2(a,t,0,f) = f(t_1,\textrm{interpolate}_1(a,t,0,f),\textrm{interpolate}_1(a,t,2,f))
\textrm{interpolate}_1(a,t,0,f) = f(t_0,a_0,a_1)
\textrm{interpolate}_1(a,t,2,f) = f(t_0,a_2,a_3)

So indeed, we get back our 2D interpolation function. The same holds for 3D.

Interpolation inside a grid

Now, say you’re in an n-dimensional grid G: \mathbb{N}^n \to V where you have values of type V at each edge G_{i_1,\ldots,i_{n}}. You’re given a point p \in \mathbb{R}^n which is inside this grid – but not neccessarily at the discrete edges – and you want to calculate the value at p using some basis interpolation function f. The first thing you need to do is calculate the grid point which is at the “lower left” of the point p. This is easy:

\tilde{p} = (\lfloor p_0 \rfloor,\ldots,\lfloor p_{n-1} \rfloor)

meaning you just round down component-wise. Next, you need all of p‘s other neighbors, and you need them in a predictable order, so our function \textrm{interpolate}_n can work on it. This sequence of neighbors is shown in figures 1 to 3 already. We have to generate the binary vector sequence of dimension n. In dimension 2, it consists of four elements:


In dimension 3, it consists of 8 elements:


Its definition is actually pretty simple: Take the numbers from 0 to 2^{n-1} and assign it the coefficient sequence (a_0,\ldots,a_n) of the 2-adic expansion (think of the binary representation of a number):

\sum_{i=0}^{n-1} a_i 2^i

Or, algorithmically, construct the sequence via (pseudocode):

vector_sequence binary_vectors(int n,vector v)
  vector_sequence result;
  if (n == 0)
    v[0] = 0;
    v[0] = 1;
    v[n] = 0;
    result += binary_vectors(n-1,v);
    v[n] = 1;
    result += binary_vectors(n-1,v);
  return result;

To calculate the interpolated point p, all you have to do is pass p - \tilde{p} as well as the binary vectors to the function interpolate and you’re done:

value interpolate_in_grid(grid g,real_vector p)
  vector ptilde = floor(p);
  // Pass the vector (0,0,...) as the starting value for binary_vectors
  vector_sequence neighbors = binary_vectors(n,null_vector<n>());

  // Add the current position to the binary vectors to get the
  // "real" neighbors of our point.
  foreach (vector &v,neighbors)
    v = ptilde + v;

  // Map from the grid _positions_ to the _values_ at the positions
  // (which could be colors, scalars, ...)
  value_sequence values;
  foreach (const vector v,neighbors)
    // Access grid at position v

      p - ptilde,
      // Start at index 0
      // Use linear interpolation
This entry was posted in Uncategorized and tagged , , , . Bookmark the permalink.

2 Responses to n-dimensional interpolation

  1. Pingback: An Introduction to Spline Modeling : Cars Blog | Everything You should Know about Cars

  2. Christ says:

    Outstanding post! I’ve been looking for a simple explanation of this for years!!! Well done!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s