Interpolation in any dimension with barycentric coordinates
Introduction
Interpolation is a method used to create a new point in a data set without doing new measurements, using the existing points. It can be used to create a continuous function from a set of points, and can have applications in various domains: it can be used to scale up the framerate of a video, to predict phenomenons with less measurement, etc. In this article, we will try to find an efficient method that can handle a lot of parameters.
Here is our problem : a company wants to determine how efficient a machine is, depending on various conditions (such as temperature, material properties, etc). The company takes multiple measurements under different conditions and then interpolates the results to guess what would happen in cases that weren't tested. The goal is not to correct the measurements or to reduce noise in the data, but simply to create a plausible continuous function from the set of points. More precisely, we will try to solve a more specific problem that is hard to solve with other methods:
The measurements can be impacted by a potentially large number of independent parameters.
The measurements aren't equally distributed: they aren't on a "grid" of predefined values, and there can be more points around some values than others.
The partial derivatives of the output aren't available.
Let's say we have n parameters. We will represent each measurement as a point in an n-dimensional space: the coordinates of the point will correspond to the value of each parameter during the measurement, and each point will be associated with its measured value. We can now try to interpolate at points where we don't have any measurements.
Visualization
To visualize how each method produces results, we will use interactive figures like the one below. It represents the problem using two parameters. Each point you see (which you can also drag around) represents a measurement. Its X and Y coordinates represent the values of its two parameters, and its color represents the measured value. Each pixel in the background indicates the result of the interpolation at its position. For this example, we used a very simple interpolation method: the output value is the value of the nearest measurement.
Nearest neighbor interpolation
Inverse distance
We can start by trying a very simple method: inverse distance weighting. If
A
0
...
A
n
are the measurement points and
α
0
...
α
n
are their values, and we want to get a result on the point M, we calculate:
Σ
i=0
n
a
i
A
i
M
Σ
j=0
n
1
A
j
M
This is a weighted average, where the weights are the inverse distances between the points. This way, the farther a point is from M, the less it will impact the result. What if a distance is zero? Then it means that M is directly on a measured point, so you can simply take its value instead.
Inverse distance interpolation
Problems with inverse distance
Although this method initially seems to work well, it has some drawbacks.
First, if there are many points in the same location (for instance, if you want more precision around a particular point), their values will disproportionately impact the result.
Bad interpolation example
On the right, four additional measurements have been made around the red point, and they impact the results across the entire space! (The color of the borders is not the same.) The impact should instead be limited to the area around the red point.
This method also has another drawback: it takes all the points into account, even when it shouldn't. In the next figure, you would expect the center of the screen to be completely blue, but the red point has impacted the result, and the center is now purplish.
Another bad interpolation example
A simpler problem
Before trying another method, we will reduce the number of points. We will first interpolate within the simplest shape possible: a triangle. If we are in a 3-dimensional space (with 3 parameters), it will be a tetrahedron, and in higher dimensions, the equivalent object. Generally speaking, in any dimension, this object is called a simplex.
We will consider a non-degenerate simplex (not completely flat, meaning the points are not all contained within the same plane, so we can't find any n−1 dimensional subspace where they all fit) in dimension n defined by the point)
A
0
..
A
n
, and try to interpolate inside it.
Fortunately, there is a nice tool for that: barycentric coordinates. Definition. (Barycentric coordinates) Given a point M,
λ
0
...
λ
n
are the barycentric coordinates of M for
A
0
...
A
n
if:
Σ
k=0
n
λ
k
A
k
M
=0
and
Σ
k=0
n
λ
k
=1 But what does it means? The first condition implies that if you consider
λ
0
...
λ
n
as weights on the points
A
0
...
A
n
, M is the barycenter of these points. The second condition is included just for convenience.
Reformulation
If we consider another point O as an origin then we have:
OM
=
Σ
k=0
n
λ
k
O
A
k
M is the weighted average of the points, and the weights are its barycentric coordinates. let el = document.querySelector("hidden#bary-reformulate"); let top = document.querySelector("hidden#bary-reformulate > .top"); top.addEventListener("click", ev => { el.classList.toggle("expanded"); }); "" For the mathematical proofs up ahead, we will provide simple explanations with interactive visualizations, as well as more formal proofs.
Formal framework
For the formal proofs, we will represent our points in an Euclidean space, meaning an affine space E, with the associated real vectorial space
E
of dimension n, containing an inner product x.y=
Σ
k=1
n
x
k
y
k
where
x
1
...
x
n
and
y
1
..
y
n
are the coordinates of the corresponding points. let el = document.querySelector("hidden#formal-framework"); let top = document.querySelector("hidden#formal-framework > .top"); top.addEventListener("click", ev => { el.classList.toggle("expanded"); }); ""
Why barycentric coordinates?
We can first make a few remarks:
If M is on
A
k
, then its coordinate for k is 1, the others are 0. For interpolation, if we use its coordinates as weights for the measurements, we will obtain the value of
A
k
, which is what we are looking for.
If M is the center of gravity of the simplex (the average of all the points), its barycentric coordinates will all be equal to
1
n1
where n is the dimension. This is also consistent.
These properties make them good candidates for interpolation. We can also prove that for all points in a simplex (that is not degenerate), these coordinates exist and are unique. We will also see later that they are easy to compute.
But first, let's introduce a basis for the space. To define a point in space, we often use an orthonormal basis, which consists of a point O as the origin of the space, and n unit vectors (where n is the dimension of the space) that are all orthogonal to each other. However, the basis does not need to be orthogonal to uniquely define each point in the space. The vectors can be of any length and point in any direction, as long as they do not all lie in the same plane.
With that in mind, we can now construct a convenient basis to define our points. First, we will take a random point of the simplex, which we will name
A
0
. It will serve as the origin. We will then name all the other points
A
1
to
A
n
.
The vectors (
A
O
A
1
,...,
A
O
A
n
) will define the basis. Because the simplex is not degenerated, we know that this basis is valid and that the vectors do not all lie in the same plane. From now on, we will refer to this basis as the "simplex basis."
We will now use the simplex basis to build the barycentric coordinates of a point M located in the simplex for each vertex. First, we know that there exist coefficients
x
k
k∈
1,n
so that:
A
0
M
=
Σ
k=1
n
x
k
A
0
A
k
(Those coefficients are the coordinates of M in the simplex basis.) Here all the vectors are relative to
A
0
, but we want them to be relative to M (to match the definition of barycentric coordinates). Thus we can add
M
A
0
on both side of the equation: (using the fact that
M
A
0
=
Σ
k=1
n
x
k
M
A
0
1
Σ
k=1
n
x
k
M
A
0
):
M
A
0
A
0
M
=
M
A
0
Σ
k=1
n
x
k
A
0
A
k
0
=
M
A
0
1
Σ
k=1
n
x
k
Σ
k=1
n
x
k
Σ
k=1
n
x
k
A
0
A
k
0
=
M
A
0
1
Σ
k=1
n
x
k
Σ
k=1
n
x
k
M
A
0
A
0
A
k
0
=
M
A
0
1
Σ
k=1
n
x
k
Σ
k=1
n
x
k
M
A
k
We can see that this resembles the definition of barycentric coordinates. We just need to set, for all k∈
0;n
y
k
=
1
Σ
k=1
n
x
k
if
k=0
x
k
otherwise
This gives us what we want :
Σ
k=1
n
y
k
M
A
k
=
0
The second condition is also verified :
Σ
k=0
n
y
k
=1 We have proven the existence of barycentric coordinates in a simplex! We have also shown that for all vertices
A
k
(except the origin), the coordinate for that vertex is equal to the coefficient
x
k
. This property will be useful later.
To prove that these coordinates are unique, here is a boring proof:
Uniqueness proof
Let M∈E, we suppose that
a
k
k∈
0,n
et
b
k
k∈
0,n
are barycentric coordinates for M in
A
k
k∈
0,n
. Then we have:
Σ
k=0
n
α
k
M
A
k
=
Σ
k=0
n
β
k
M
A
k
This implies:
Σ
k=0
n
M
A
k
α
k
β
k
=0
Σ
k=0
n
A
0
M
α
k
β
k
Σ
k=0
n
M
A
k
α
k
β
k
=0
Σ
k=0
n
A
0
A
k
α
k
β
k
=0 But
A
0
A
1
...
A
0
A
n
form a basis of
E
, so for all k∈
0;n
,
α
k
β
k
=0 which means
a
k
=
b
k
let el = document.querySelector("hidden#existence-proof2"); let top = document.querySelector("hidden#existence-proof2 > .top"); top.addEventListener("click", ev => { el.classList.toggle("expanded"); }); ""
Barycentric coordinates in a simplex
We will now try to find a simple (and easy to compute) expression for barycentric coordinates. To make things simpler, we will assume that M is inside the simplex. As we saw before, the barycentric coordinate of a vertex
A
k
is equal to the coefficient
x
k
of the vector to this vertex
A
0
A
k
in the simplex basis. Thus, we just have to find an expression for that coefficient. If we take all the other vertices of the simplex, we know that they define a plane
P
k
, and that
A
k
is the only point out of it. This means that
A
0
A
k
is the only vector not perpendicular to
n
, the unit normal vector of
P
k
pointing towards
A
k
. We have:
A
0
M
.
n
=
Σ
i=1
n
x
i
A
0
A
i
.
n
=
x
k
A
0
A
k
.
n
=
x
k
A
0
A
k
.
n
A
0
A
k
.
n
is just the height of the simplex and
A
0
M
.
n
is just the height of M. Thus we have:
x
k
=
Height of
M
Height of
A
k
Just below is a visualization of the situation for a 2-simplex (a triangle). The blue lines are the vector representing M in the simplex basis, and the green line represents the heights of M, the yellow line the height of
A
k
.
x
k
=
Height of
M
Height of
A
k
=
999
999
=999
Coefficient visualisation
Same proof, but more rigorous
Let set for all k∈
1,n
,
b
k
=
A
0
A
k
.
b
1
...
b
n
\
b
k
defines a hyperplane
B
k
of
E
. We will name
P
k
the associated plane in E. Let H be the orthogonal projection of M on
B
k
and K the orthogonal projection of
A
k
on
B
k
.
Let k∈
1,n
. We have :
A
0
M
=
Σ
i=1
n
x
i
A
0
A
i
A
0
M
=
Σ
i∈
1,n
\
k
n
x
i
A
0
A
i
x
k
A
0
A
k
Where
x
i
i∈
0;n
are the barycentric coordinates of M in
A
0
...
A
n
Let
n
the normal vector of
B
k
pointing towards
A
k
. We have:
x
k
A
0
A
k
.
n
=
HM
x
k
A
0
A
k
.
n
=
HM
x
k
K
A
k
.
n
A
0
K
.
n
=
HM
Since K and
A
0
are both in
P
k
, and
K
A
k
and
HM
are both normal to
B
k
.
x
k
K
A
k
=
HM
This gives us the expression of
x
k
:
x
k
=
HM
K
A
k
Now by following the proof of the existence of the coordinates, we have
λ
k
=
x
k
. let el = document.querySelector("hidden#expression-proof2"); let top = document.querySelector("hidden#expression-proof2 > .top"); top.addEventListener("click", ev => { el.classList.toggle("expanded"); }); ""
Area and volume
To better visualize how barycentric coordinates behave, we can use another formula. In 2 dimensions, the area of a triangle is given
basis
height
2
. We can reformulate the expression: Let l≠m≠k
x
k
=
Height of
M
Height of
A
k
=
A
A
l
A
m
M
A
A
l
A
m
A
k
This is the ratio of the triangle's area where
A
k
have been replaced by M, by the area of the original triangle. Here is an example:
blue area
gray area
=
999
999
=999
Barycentric coordinates with area
The same concept applies in three dimensions: the coordinates are the ratio of the volumes.
blue volume
total volume
=
999
999
=999
Barycentric coordinates with volume
It makes sense because the sum of the small areas/volumes will be the total simplex area/volume, so the sum of the coordinates will be 1.
However, this formula is only useful for visualization, since heights are easier to compute. The program used to compute barycentric coordinates in the next interactive figures will only use the heights.
You can see that this method works well for interpolation:
Barycentric interpolation in one simplex
Since we have maintained a general approach throughout, this method works in any dimension!
Interpolation in a cloud of points
Even tho barycentric interpolation works very well, the case presented above isn't very practical : we only have n points, and the area where we can interpolate is very small. Ideally, we would like to interpolate within any cloud of points over a larger area.
We can try to triangulate the space with those points (by choosing simplexes that cover the entier space), but there are some drawbacks:
First, there is no unique way to achieve this (see figures below).
Triangulation algorithms are already complex in two dimensions; extending them to higher dimensions is likely to be very challenging.
If interpolation is needed at a single point, there is no need to compute a triangulation for the entire space—it would be more efficient to consider only the simplexes surrounding the point.
Two different interpolations
Additionally, bad triangulations can lead to terrible interpolation results.
A terrible triangulation
Another approach is to select the n nearest points to form a simplex. However, there are cases where the formed simplex does not contain the point of interest.
The white point is not in the triangle formed by the 3 nearest points
The final method
To get better results, we can try interpolating in multiple simplexes. We can select an arbitrary number of simplexes that contains the point, and calculate the average of their interpolations. To choose the best simplexes, several criterias can be considered:
The smallest volumes
The smallest surfaces
The smallest
surface
volume
ratio (to avoid stretched simplexes).
The choice of criteria and the number of simplexes to consider should be adapted to the nature of the interpolation problem: some datasets can be more adapted to a specific method.
Here is an example with the volume criterion, with 2 simplexes used to interpolate:
Average of interpolations of the 2 smallest simplexes
As you can see, the algorithm creates discontinuities around points that are near the border of a triangle. Also, the algorithm seems to try interpolating in stretched triangles, resulting in the blue stripes and purple stripes around the center.
To address this issue, we can use a weighted average. Let S be the set of all simplexes that contains the point to interpolate. Define w a function that takes a simplex and returns a criterion as described above, and i a function that returns the interpolated value in the simplex. We can then compute:
Σ
S∈A
i
s
w
s
n
Σ
S∈A
1
w
s
n
Here, w returns the surface of the simplex, and n=2:
Interpolation with average of multiple simplexes
The result may not be perfect, but if discontinuities are a significant issue, it is possible to smooth the results further by calculating the interpolation at multiple points around the desired value.
Conclusion
Even though we focused primarily on 2D examples (particularly to demonstrate the value of interpolation points with colors), one of the greatest advantages of this method is its ability to be generalized to any dimension. We choose a very specific problem at the beginning, but the way we used barycentric coordinates in simplexes around the point we interpolate makes this method very flexible compared to the other methods we introduced earlier. It allows you to decide how to aggregate interpolation values from different simplexes to obtain the final result, as demonstrated in the last section.
There are still some drawbacks to this method. It is computationally more intensive than the simpler methods we discussed earlier (although it remains reasonable) and can only interpolate within the convex hull of the points (i.e., the interior of the shape formed by connecting all the outer points of the cloud).
Links and source code
If you are really interested, here is an implementation of this method, as well as the source code of the website: Github repository
This article presents a method for obtaining barycentric coordinates in shapes more complex than simplexes, which could be a great solution for improved results: a great solution to get better results: Error Estimates for Generalized Barycentric Interpolation