Example of how to interpolate a Z value for a point using the Delaunay Triangulation in JTS

Question

This is a fairly remedial question. I have been looking at the documentation for the JTS DelaunayTriangulationBuilder and I am at a loss as to how to do what seems like it should be a simple thing. I wish to take a series of points, triangulate them, and then interpolate the Z value of a random point within that mesh. It's non-obvious from a cursory reading how to do this. Any ideas?

Answer 1

After you've loaded up the triangulation object, call getSubdivision() on it to get the triangulation. It uses a quad-edge data structure, which you'll need later. (It's easier to understand if you know what a half-edge or winged-edge representation is.) The resulting QuadEdgeSubdivision has a method locate that, given a coordinate, returns one of the edges of the enclosing triangle (as a quad-edge). Get its origin vertex with orig() and its destination vertex with dest(). Get another edge with oNext() Its destination vertex is the third vertex (also dPrev().origin() is the same vertex). Now that you have the three vertices, represent your test point as a Vertex and call interpolateZValue.

For example:

public static double 
interpolateZ(DelaunayTriangulationBuilder triangulation,
             Coordinate coordinate) {
    QuadEdgeSubdivision quadEdgeSubdivision = triangulation.getSubdivision();
    QuadEdge edge = quadEdgeSubdivision.locate(coordinate);
    return new Vertex(coordinate.x, coordinate.y)
            .interpolateZValue(edge.orig(), edge.dest(), edge.oNext().dest());
}

You're right, though. It's not obvious how to do this from reading their API.

Answer 2

I'm not familiar with JTS DelauneyTriangulationBuilder, but it sounds like you have a collection of points (x,y,z), and you are submitting the 2D (x,y) pairs to the triangulator. This gives you a planar triangulation of the (x,y) points, but also a mesh who's vertices are the original (x,y,z) points.

Once you have a triangulation, you wish to find the point (p,q,r) on the mesh that corresponds to the planar point (p,q). To do so, find the triangle T of the Delauney triangulation that contains (p,q). Find the barycentric coordinates of (p,q) relative to T, and use these to compute a weighted average r of the z values corresponding to the vertices of T. That weighted average is the Z value you're looking for. In other words, (p,q,r) is on the mesh.