Bilinear interpolation on integer coordinates within a Delaunay triangulation

Question

I have a planar Delaunay triangulation consisting of about 1 million triangles. Each vertex is tagged with several scalar metrics [1], and I would like to see a fast, simple interpolation of each of those metrics on the same regular grid. For reference, the union of my triangles covers about 10 million grid cells having (integer) coordinates. [2]

When I say simple, I mean simple. Bilinear is fine! My understanding is that this is (a) basically what GPUs do for a living, and (b) probably the subject of innumerable homework assignments. I myself am a government researcher in public health, so it's not homework for me. :-)

In my slow but correct reference implementation, I can compute the following in about 10 minutes:

For each triangle T:

1. The set G of all (integer) Cartesian coordinates within the bounding box of T;
2. Barycentric coordinates (u, v, w) for each (x, y) in G;
3. Rejection of (u, v, w) that are not all positive — that is, inside T;
4. The weighted sum (uz_1 + vz_2 + w*z_3) for each remaining coord in T, where z_1, z_2, and z_3 are, for a given metric [1], the scalar values at the vertices of T.

I really need steps 1-3 to be fast; step 4 is trivial but it's my end goal. Ideally the solution would take either of the following forms:

• A suitably licensed (GPL is OK) library with a dead simple API; or
• An explanation that's clear enough that it's obvious how an intermediate programmer could code it up in Fortran, R, Python, or C.

A classic formulation of this task is the "TIN to DEM" terrain-modeling job. But it seems the reverse is more commonly needed these days (?)

Some basic cleanup, like removing duplicates when a point falls exactly on an edge or vertex shared by 2+ triangles, is OK too.

Many thanks in advance for your time and attention. I'll clean up formatting and edit per suggestions once I'm off the train!

Footnotes:

[1] Elevation, temperature, and humidity. [2] Integers in the sense that they're spaced 20x20m apart on a UTM grid. So just scale by 20.

Answer 1

Though I don't see anything that might explain the slowness in your description, here is how I would handle it. The basic ingredient is indeed the "triangle scanner".

Start by sorting the three vertices on Y. This takes three comparisons and there are just six possible configurations. Loop on the integer ordinates from top Y to middle Y then middle Y to bottom Y. For every ordinate, the intersections with the left and right sides gives you an interval.

Loop on the integer abscissas in that interval, from left to right. The double loop will only visit grid nodes that belong to the triangle.

Rather than using barycentric coordinates, you can establish the equation of the plane, i.e. evaluate the coefficients of Z = a X + b Y + c, and use this formula for interpolation. (You can even compute the values incrementally, i.e. Z(X + 1) = Z(X) + a, but for a small number of points per triangle, I am not sure this is worth.)

It is easy to avoid duplicate points by relying on a simple convention: only produce points that fall on the left side of the triangle, not those that fall on the right side (these will be produced by the triangle to the right.)

Some car must be exercised to handle special cases, such as an horizontal side at integer ordinate, but this is manageable.

The total workload will be sensitive to the number of triangles, the Y extent of the domain and the number of covered grid nodes, counting a handful of arithmetic operations for each of these factors. For one million triangles and ten million grid cells, running times below a second are not unlikely.