Linear interpolation given a triangulation in arbitrary dimensions

Specific Domains Numerics
1

Hi all,

I’m in need of the equivalent of LinearNDInterpolator in scipy for julia, in order to do interpolation of irregularly sampled points (not on a grid). There are existing packages such as NaturalNeighbours.jl, but any implementation I’ve seen only works on 2d data so was trying to do my own. Posting here in case someone already has experience with something similar and can share some tips (or perhaps point me to an existing implementation I’ve missed).

Idea is to do a Delaunay triangulation in arbitrary dimensions and then apply barycentric interpolation, same as the scipy function linked above. To do this the complex parts are:

  • Constructing the Delaunay triangulation in arbitrary dimensions
  • For a given point efficiently finding the vertices of the “triangle” that contains it.

After that the barycentric interpolation is a straightforward step. The first part can be quickly achieved with the Quickhull.jl package:

using Quickhull
npoints = 1000
points = randn(3,npoints)
tri = delaunay(points)
facets(tri)

But from this point forward I end up being stuck in terms of existing packages. I could go ahead and try to implement the rest, but I feel I might just be missing existing (and efficiently implemented) ways to achieve the remaining step of finding the vertices for the triangle surrounding a given point.

Thanks in advance for any help!

2

Do you really need a triangulation or do you only need the interpolation in the end? If you want to interpolate irregularly sampled points in an arbitrary dimension, I would recommend using a method that’s well suited for scattered data interpolation. One prominent example are radial basis functions. I have the package KernelInterpolation.jl, which provides an implementation of scattered data interpolation in arbitrary dimension. If you need a linear interpolation, you can try kernel = PolyharmonicSplineKernel{NDims}(1), where NDims is there dimension of the samples. However, if you only need an interpolation, not necessarily a linear one, I would recommend another kernel, such as ThinPlateSplineKernel.

2 Likes
3

Hi Joshua,

thanks for pointing out your package, had not seen that one. I do not really need a triangulation, its more that the requirements I have are that:

  • the interpolation is continous
  • it recovers the input values at the data points
  • speed of each interpolation must be very fast (a possibly extensive initialization does not matter). Runtime of interpolations should not scale steeply with the number of data points used.

All of which can be satisfied with linear interpolation based on a triangulation. I’m a bit unfamiliar with how kernel interpolation works, but can it satisfy those requirements?

4

I was working on a similar problem. My code is in a private repo for now and still a bit messy, but I’d be happy to share it if you’re interested.

For a given point efficiently finding the vertices of the “triangle” that contains it.

I didn’t need this to be particularly efficient, since I was working on estimation problem where I repeatedly evaluate changing interpolators on the same data points, but I simply used the following to check whether a point is in a “triangle”

struct Simplex{P,T}
  vertices::P
  M::T
end

function Simplex(points)
  Simplex(points, inv(vcat(ones(eltype(points),size(points,2))', points)))
end


function barycentric(point, triangle::Simplex)
  @views λ = triangle.M[:,1] + triangle.M[:,2:end]*point
  return (λ)
end

function isin(point, triangle::Simplex)
  @inbounds for i in 1:size(triangle.M,1)
    λ = triangle.M[i,1]
    @inbounds for j in 2:size(triangle.M,2)
      λ += triangle.M[i,j]*point[j-1]
    end
    if (λ<-eps(eltype(point))) #0)
      return(false)
    end
  end
  return(true)
end

and looped over all triangles. A faster way that I considered but didn’t implement would be to keep track of the location of the triangles relative to each other so that a positive or negative barycentric coordinate for one triangle eliminates the need to check a bunch of others.

5

The first two points are satisfied for the kernel based interpolation. In fact, the smoothness of the resulting interpolant depends on the smoothness of the kernel function, but all of them are at least continuous. Regarding the performance, I cannot guarantee it’s the best as it needs to solve a linear system of the size of the number of input nodes, but maybe you could give it a try? It should be quite straightforward once you have the nodes and function values at the nodes. If it is not fast enough, one could try smarter ways to solve the linear system I would consider to implement if required.

6

Following @Paul_Schrimpf 's example, ended up making an implementation that just brute forces through the simplexes. This way I could compare performance against KernelInterpolations. Made a gist with it here, and it allows me to do this:

# create test data
test_function = (x) -> sin(x[1]^2 + x[2]^2)
npoints = 1000
points = rand(2,npoints)
values = zeros(npoints)
for i in 1:npoints
    values[i] = test_function(points[:,i])
end
# create interpolant
si = SimplexInterpolant(points);

# evaluate interpolant in regular grid
xs = range(0, 1, length=100)
ys = range(0, 1, length=100)
zs = [0.0 for x in xs, y in ys]
for (i,x) in enumerate(xs)
    for (j,y) in enumerate(ys)
        info = interpolation_info([x,y], si)
        zs[i,j] = compute_simplex_interpolation(info, values)
    end
end

#plot result
using CairoMakie, GeometryBasics
fig, ax, hm = heatmap(xs, ys, zs)
tri = delaunay(points)
wireframe!(ax, GeometryBasics.Mesh(tri))
Colorbar(fig[:, end+1], hm)

save("plot.png", fig)

fig

plot
plot1200×900 526 KB

So with this example I did some benchmarking against KernelInterpolations. Using the same data initialized above I have:

using KernelInterpolation

nodeset = NodeSet(transpose(points))
kernel = PolyharmonicSplineKernel{dim(nodeset)}(1)
ff = test_function.(nodeset)
itp = interpolate(nodeset, ff, kernel)
##
using BenchmarkTools
@benchmark itp(pt) setup=(pt=rand(2))

##
@benchmark begin
    info = interpolation_info(pt, $si)
    compute_simplex_interpolation(info, $values)
end setup=(pt=rand(2))

which gives me a median runtime of 41 microseconds for KernelInterpolations and 13 microseconds for my interpolation. Both seem to scale more or less linearly with number of points which I pretty much expect in my simplex interpolation. But an enormous amount of the time I spend is just in finding the appropriate simplex, would not be surprised if I could get a ~2 order of magnitude speedup just by implementing an efficient algorithm to find the right simplex.

7

Oh, I see. You care about evaluating the interpolant being fast, not computing the interpolation? I previously thought of the performance of the interpolation process, which is the computationally expensive part.
I just checked with ProfileView.jl, where most of the time is spent in evaluating the linear combination in KernelInterpolation.jl and all the time is basically only spent in computing differences x - y of points and computing norms. I found that by using some simple tricks I could reduce a median of 57 microseconds on my machine for your example down to 16 (see PR). So performance seems to be similar in this case compared to your implementation. The new version (v0.3.1) should be available in a couple of minutes.

8

Yes, just a kd tree should work wonders here. However, most of the kd trees I’m aware of in Julia right now search points, whereas you want to search simplices, and the two problems are not quite equivalent (a simplex can lie on both sides of a separating plane). Fortunately, implementing this is not too hard, and it would be useful to have for finite-element packages — I linked some sample code in: Error with the interpolation of 3D model · Issue #831 · gridap/Gridap.jl · GitHub (and maybe one of the many FEM packages has already done it too?).

9

Ah interesting! I hadn’t thought of building a kd tree of bounding boxes / geometries before for that reason, but that could also be something to do.

In geospatial we encounter the problem of locating which geometry a point is in quite a lot. To that end there are packages like SpatialIndexing.jl (should work in n dimensions with bounding boxes) and SortTileRecursiveTree.jl (only works in 2d for now, but could be generalized) that can help you do this. The way you set this up is you construct the tree with the bounding boxes, and then query the tree to get the potentially intersecting bounding boxes. After that you can brute force through the much smaller collection of simplices.

An interesting algorithm at scale which is enabled by GeometryOps.jl and its SpatialTreeInterface (but you can do this with AbstractTrees.jl and SortTileRecursiveTree.jl too) is to consider your grid to be a “tree” of its own, if it is a regular (or even irregular) grid. Then, you can iterate down both trees at once in order to decrease the amount of time you spend searching.

If you have a strong usecase for that then I can help out with a PR to make SortTileRecursiveTree.jl n-dimensional.

1 Like