Monday, 1 November 2021

JTS Polygon Triangulation, at last

A (long) while ago I posted about "soon-to-be-released" JTS code for polygon triangulation using Ear Clipping.  It turned out it was actually in the category of "never-to-be-released".  However, later  I worked with a student, Dan Tong, on a coding exercise sponsored by Facebook. We decided to tackle polygon triangulation.  By the end of the project he implemented a functional algorithm, including handling holes (via the hole-joining technique described by Eberly).  To make it release-ready the code needed performance and structural improvements.  This has finally happened and the code is out!  

Polygon Triangulation

The new codebase has been substantially rewritten to improve modularity and the API.  A new data structure to represent triangulations has been introduced, called a Tri  It is simpler, more memory efficient, and easier to work with than the QuadEdge data structure previously used in JTS. It provides an excellent basis for developing further algorithms based on triangulations (of which more below)

Most of the reworking involved implementing key performance optimizations, utilizing existing and new JTS spatial index structures: 

  • the Hole Joining phase now uses a spatial index to optimize the test for whether a join line is internal to a polygon
  • the Ear-Clipping phase has to check whether a candidate ear contains any vertices of the polygon shell being clipped.  This has been made performant by using a new spatial index structure called the VertexSequencePackedRtree

Performance must be benchmarked against an existing implementation.  Geometry maestro Vladimir Agafonkin has developed an Ear Clipping algorithm in Javascript called Earcut.  It's claimed to be the fastest implementation in Javascript.  Indeed the performance comparison shows Earcut to be over 4 times faster than some other Javascript triangulation implementations.  Earcut has been ported to Java as earcut4j, which allows comparing JTS to it.   

DatasetSizeJTS TimeEarcut4j Time
Lake Superior3,478 vertices / 28 holes18 ms13 ms
Canada10,522 vertices35 ms26 ms
Complex polygon44,051 vertices / 125 holes1.1 s0.64 s
North America with Lakes115,206 vertices / 1,660 holes13.3 s5.4 s

So JTS ear clipping is not quite as fast as Earcut4j.  But it's still plenty fast enough for production use!

The triangulation algorithm is exposed as the PolygonTriangulator class (deliberately named after the effect, rather than the algorithm producing it).  Here's an example of the output:

Ear-Clipping Triangulation of Lake Superior

Constrained Delaunay Triangulation

The output of PolygonTriangulator is guaranteed to be a valid triangulation, but in the interests of performance it does not attempt to produce a high-quality one.  Triangulation quality is measured in terms of minimizing the number of "skinny" triangles (or equivalently, maximizing the sum of the interior angles).  Optimal quality is provided by Constrained Delaunay Triangulation (CDT).  In addition to being optimal in terms of triangle shape, the CDT has the nice properties of being (essentially) unique, and of accurately representing the structure of the input polygon.  It also removes the influence that Hole Joining holes has on the raw triangulation provided by Ear-Clipping.

I originally proposed that a CDT can be achieved by a Delaunay Improvement algorithm based on iterated triangle flipping.  It turned out that this technique is not only effective, but also pleasingly performant.  This is implemented in the ConstrainedDelaunayTriangulator class.  Running it on the example above shows the improvement obtained.  Note how there are fewer narrow triangles, and their arrangement more closely represents the structure of the polygon.

Constrained Delaunay Triangulation of Lake Superior

This code will appear in the next JTS version (1.18.3 or 1.19), and has already been released in GEOS 3.10.

Future Work

It seems possible to adapt the machinery for Ear Clipping to provide two other useful constructions:

  • Polygon Concave Hull ("Outer Hull") - this requires a triangulation of the area between the polygon and its convex hull.  That can be obtained by running Ear Clipping outwards, rather than inwards. A Concave Hull can be constructed by limiting the number of triangles generated based on some criteria such as area, longest edge, or number of vertices.  The result is a simplified polygon with fewer vertices which is guaranteed to contain the original polygon. 
  • Polygon Inner Hull - the triangulation constructed by Ear Clipping can be "eroded" to provide a simpler polygon which is guaranteed to be contained in the original.
The property of the Constrained Delaunay Triangulation of representing the structure of the input polygon could provide a basis for computing the following constructions: 
  • Approximate Skeleton construction (also known as the Medial Axis) by creating edges through the CDT triangles
  • Equal-Area Polygon Subdivision, by partitioning the graph induced by the CDT
The Tri data structure provides a simple representation of a triangulation.  I expect that it will facilitate the development of further JTS triangulation algorithms, such as:
  • Concave Hull of Point Set
  • Constrained Delaunay Triangulation of Lines and Points