0

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.

dholstius
  • 1,007
  • 1
  • 7
  • 16
  • Quick questions: What precision do you need for the output and, if it's not excessive, can you just use, say, OpenGL/OpenGL-ES on a GPU? As you say, this is exactly the sort of thing GPU's do (and do exceedingly rapidly). – Simon F May 12 '16 at 15:27
  • I only need precision of about 1% (maybe 0.1%) of the full scale of each metric. So for example, if the metric "z" were temperatures, and they ranged from 0 to 100 C over the full domain, then I only need interpolations of z to a precision of about 1.0 or maybe 0.1 C. Many triangles will be "flat", having vertices that are all very close to the same z. I've tried to generate the triangulation such that the gradient will be reasonably smooth — that is, where I expect the gradient to be steep, I've generated a finer mesh. I could plot a histogram of the range of z by size of triangle ... ? – dholstius May 12 '16 at 15:34
  • This looks promising, except that I already have the Delaunay triangulation that I want to use: http://rncarpio.github.io/delaunay_linterp/ – dholstius May 12 '16 at 15:47
  • Also promising: https://github.com/geotrellis – dholstius May 12 '16 at 16:43
  • Looks like what I want: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html – dholstius May 12 '16 at 17:46
  • Ooh: http://www.r-bloggers.com/barycentric-interpolation-fast-interpolation-on-arbitrary-grids/ – dholstius May 12 '16 at 17:50
  • Ten minutes sounds really slow. I would expect a C compiled implementation to run in a few seconds at most. What language/environment are you using ? –  May 12 '16 at 19:56
  • I'm using R. Might be a memory I/O issue (my naive implementation uses around 10 GB, which a good implementation shouldn't, since I keep copies of intermediate steps). Might be overhead from the interpreter. Etc. I used the "profvis" package to eliminate some obvious bottlenecks, and that's how I got it down to 10 minutes. – dholstius May 12 '16 at 21:55
  • @dholstius: with more care, you can probably stay below 1 GB. –  May 12 '16 at 22:03

1 Answers1

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.)

enter image description here

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.