## 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

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

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

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

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})$

with

$\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

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:

$((0,0),(0,1),(1,0),(1,1))$

In dimension 3, it consists of 8 elements:

$((0,0,0),(0,0,1),(0,1,0),(0,1,1),(1,0,0),(1,0,1),(1,1,0),(1,1,1))$

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;
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 $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
values.push_back(g[v]);

return
interpolate(
n,
values,
p - ptilde,
// Start at index 0
0,
// Use linear interpolation
linear_interpolation);
}

This entry was posted in Uncategorized and tagged , , , . Bookmark the permalink.

### 2 Responses to n-dimensional interpolation

1. Christ says:

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