2D irregular interpolation

I have worked out a scheme for linear interpolation of functions `f(x,y)` that I think might be fast enough to be useful (and might be possible to optimize really well with the same metaprogramming tricks we use on B-splines). It's probably not novel, but I don't know what it's called, so I'm just going to describe it.

Given lists of points `(x,y,z)` (maybe on the form of equal-length lists `xs`, `ys` and `zs`) representing a discretization of a function `z = f(x,y)`, we can find the value of `z` at an arbitrary point inside the domain using the following steps:
- Initialization:
- Create a triangular mesh on the xy-plane
- Evaluation:
1. Find the three corners of the triangle that contains the sought point `(x,y)`. Let's call them `P`, `Q` and `R`.
2. Evaluate `z` as follows:

```
A = cross(Q-P, R-P)
d = dot(A, P)
z = (d - (x * A[1] + y * A[2])) / A[3]
```

(These relations follow from finding the equation of the plane through the three points, and then solving for `z`.)

I would love to have this functionality in the package, and I could take a stab at an implementation (at least as a POC) but I have no idea how to build the triangular mesh. Maybe a [Delaunay triangulation](https://en.wikipedia.org/wiki/Delaunay_triangulation) is what we want? Is there a package that does that?