Barycentric interpolation generalises linear interpolation to arbitrary dimensions. It is very fast although suboptimal if the function is smooth. You might now it as algorithm 21.7.1 in Numerical Recipes (Two-dimensional Interpolation on an Irregular Grid). Using package geometry it can be implemented in a few lines of code in R.
Here’s a quick explanation of what the algorithm is doing: when interpolating to approximate the value of a function at , you want to see what the value of
is at neighbouring grid points and return some kind of weighted average of these values (with higher weights given to closer grid points). We start with random grid points:
X <- cbind(rnorm(n),rnorm(n))
Barycentric interpolation on the plane uses just three near-neighbours to interpolate. As a first step a triangulation is built: usually but not necessarily the Delaunay triangulation. This divides the convex hull of the grid into triangular tiles:
dn <- delaunayn(X)
The function delaunayn returns the index of the three points that make up each triangle. The tiles define which neighbours to use in the interpolation: for example if we interpolate anywhere within the blue tile we will be using the measured value of at the three highlighted points.
The only thing that we still need to define is the relative weights of the three points in the interpolated value. This is where barycentric coordinates come in: given a point in the interior of a triangle formed by
, we can express
as a weighted average (convex combination) of
:
where and
. The
‘s are called the barycentric coordinates of
, and we can use them as weights in the interpolation:
For a linear function this interpolation is exact, so barycentric interpolation makes sense for functions that are roughly piecewise linear on the triangulation.
The function tsearch will assign each interpolated location to its Delaunay tile, and return its barycentric coordinates. To avoid slow R loops I use a sparse matrix trick to compute the interpolation. The whole thing is essentially five lines of code:
#2D barycentric interpolation at points Xi for a function with values f measured at locations X #For N-D interpolation simply replace tsearch with tsearchn and modify the sparse matrix definition to have non-zero values in the right spots. interp.barycentric <- function(X,f,Xi) { require(geometry) require(Matrix) dn <- delaunayn(X) tri <- tsearch(X[,1],X[,2],dn,Xi[,1],Xi[,2],bary=T) #For each line in Xi, defines which points in X contribute to the interpolation active <- dn[tri$idx,] #Define the interpolation as a sparse matrix operation. Faster than using apply, probably slower than a C implementation M <- sparseMatrix(i=rep(1:nrow(Xi),each=3),j=as.numeric(t(active)),x=as.numeric(t(tri$p)),dims=c(nrow(Xi),length(f))) as.numeric(M%*%f) }
In what package is ‘delaunayn’? I wish all coders would list ALL packages that are needed to implement their sample scripts; 90% of the authors DON’T!!
In package “geometry”. You can usually find out easily by typing in the function’s name in rseek.org.
Hello! How did you plot the grid? It looks really nice.
never mind, of course as soon as i commented I found out about the tetramesh package, thanks for the post, though, very useful!
You’re welcome!
Hi,
Here is some friendly advice from a fellow user of this function. If you get the following cryptic error message, it probably means that you are trying to interpolate to points that lie outside the “convex hull” (outermost polygon that can be defined from the data). Go back and check your points Xi.
Error in sparseMatrix(i = rep(1:nrow(Xi), each = 3), j = as.numeric(t(active)), :
NA’s in (i,j) are not allowed
Now I have need to do this kind of interpolation in 3d. I see that function delaunayn() works for higher dimensions, but tsearch() from the geometry package is only for 2d (x,y). Is there a way to make this work for a 3d problem?
Nevermind, I found the comment by Simon about using tsearchn. Here is a version of the function for higher-order problems. Note the other change needed is to the “each” parameter for argument i in the sparseMatrix() call.
interp.barycentricN <- function(X, f, Xi) {
require(geometry)
require(Matrix)
dn <- delaunayn(X)
tri <- tsearchn(X, dn, Xi)
# For each line in Xi, defines which points in X contribute to the interpolation.
active <- dn[tri$idx,]
# Define the interpolation as a sparse matrix operation. Faster than using apply, probably slower than a C implementation.
M <- sparseMatrix(i=rep(1:nrow(Xi), each=ncol(X)+1), j=as.numeric(t(active)), x=as.numeric(t(tri$p)), dims=c(nrow(Xi), length(f)))
as.numeric(M%*%f)
} # end interp.barycentricN()
# test
f <- c(1,1,1,2,2,2,3,3,3,1,1,1,2,2,2,3,3,3,rep(4,9))
ijk <- expand.grid(i=1:3, j=1:3, k=1:3)
X <- as.matrix(ijk) + runif(n=27, min=-0.1, max=0.1)
Xi <- matrix(c(1.1,1.3,1.2,
1.7,1.8,2.8,
1.2, 1.7, 2
), byrow=T, ncol=3)
interp.barycentric3d(X, f, Xi)
[1] 1.365510 3.657084 1.907660
# Note I used irregular coordinates for X. If you have a regular grid, tsearchn() will give # a warning that there are Degenerate simplices.