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).

We first consider

*linear*interpolation in

*one*dimension. In this case, you have an vector space of values – meaning the elements have the following two operators:

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 and a scalar and calculate

In figure 1, this is depicted for . Now suppose you want to interpolate not on a line but on a square. In two dimensions, you have four points instead of two (the number always doubles) and two scalars (or one point ), see figure 2 for an example.

To interpolate the point between the four points, you have to do interpolations, two for each line and then “inbetween” the lines.

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.

The corresponding interpolations are:

## 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 of points and one “position” vector , define the interpolation in n dimensions as:

with

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

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.

To show that the formula makes sense, let’s apply it to the 2D case: We’ve got 4 values , a point and an interpolation function . The first function call has to be with . We get:

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 where you have values of type at each edge . You’re given a point which is inside this grid – but not neccessarily at the discrete edges – and you want to calculate the value at using some basis interpolation function . The first thing you need to do is calculate the grid point which is at the “lower left” of the point . This is easy:

meaning you just round down component-wise. Next, you need all of ‘s other neighbors, and you need them in a predictable order, so our function 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 to and assign it the coefficient sequence of the *2-adic expansion* (think of the binary representation of a number):

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

vector_sequence binary_vectors(int n,vector v) { vector_sequence result; if (n == 0) { v[0] = 0; result.insert(v); v[0] = 1; result.insert(v); } else { 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 , all you have to do is pass 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 values.push_back(g[v]); return interpolate( n, values, p - ptilde, // Start at index 0 0, // Use linear interpolation linear_interpolation); }

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