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

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)

}