Sunday, 17 May 2020

JTS Overlay - the Next Generation

In the JTS Topology SuiteOverlay is the general term used for the binary set-theoretic operations intersection, union, difference and symmetric difference.  These operations accept two geometry inputs and construct a geometry representing the operation result.  Along with spatial predicates and buffer they are the most important functions in the JTS API.

Intersection of MultiPolygons

Overlay operations are used in many kinds of spatial processes. Any system aspiring to provide full-featured geometry processing simply has to provide overlay operations.  In fact, many geometry libraries exist solely to provide implementations of overlay. Notable libraries include the ESRI Java API, Clipper and wagyu.  Some of these provide overlay only for polygons, which is the most difficult case to compute.

Overlay in JTS 

The JTS overlay algorithm supports the full OGC SFS geometry model, allowing any combination of polygons, lines and points as input.  In addition, JTS provides an explicit precision model, to allow constraining output to a desired precision.  The overlay algorithm is also available in C++ in GEOS.  There it provides overlay operations in numerous systems, including PostGIS, QGIS, Shapely, and r-sf.  This codebase has had an long lifespan; it was developed back in 2001 for the very first release of JTS, and while there have been improvements over the years the core of the design has remained unchanged.

However, there are some long-standing issues with JTS overlay.  The most serious one is that in spite of much valiant effort over the years, overlay is not fully robust.  The constructive nature of overlay operations makes them particularly susceptible to the robustness issues which are notorious in geometric algorithms using floating-point numerics.  It can happen that running an overlay operation on seemingly innocuous, valid inputs results in the dreaded TopologyException being thrown.  There is a steady trickle of issue reports about this in JTS, and even more for GEOS (such as here, here and here...).

Another issue is that the codebase is complex, and thus hard to debug and modify. Partly this is because of the diversity of inputs and the explicit precision model.  To support this the JTS overlay algorithm has a rich and detailed semantics.  But some of the complexity is due to the original design of the code. This makes it difficult to incorporate new ideas for improvements in performance and robustness.

Next Generation Overlay

So for many years it's been on my mind that JTS overlay needs a thorough overhaul. I chipped away at the problem over time, but it was clear that it was going to be a major effort.  Now, thanks to the support of my employer Crunchy Data, I've at last been able to focus on a complete rewrite of the JTS overlay module.  It's called OverlayNG.

The basic algorithm remains the same:
  1. Extract the input linework, and node it together
  2. Build a topology graph from the noded linework
  3. Compute a full topological labelling of the graph
  4. Extract the resultant polygons, lines and points from the graph
This algorithm is time-tested and is able to handle the complexities of multiple geometry types and topology collapse. The new codebase benefits from 20 years of experience to become simpler and more modular, with increased testability, and potential for reuse.

OverlayNG has the following improvements:
  • A snap-rounding noder is available to make overlay fully robust.  This eliminates the possibility of TopologyExceptions (when an appropriate precision model is used).
Intersection operation with Snap-rounding 
and Topology Collapse removal
  • Snap-rounding allows full support for specifying the output precision model.  The precision model can be specified independently for each overlay call, which is more flexible and easier to use.  The use of snap-rounding also provides fully valid precision reduction for geometries.  This makes it feasible for the first time to fully operate in a fixed-precision regime.
Precision Reduction turned all the way up to 11
  • Significant performance optimizations are included (notably, one which makes polygon intersection much faster in many cases)
Intersection of a MultiPolygon with a grid (7x faster with OverlayNG)
  • Pluggable noding allows providing different noding strategies.  One use is to run OverlayNG with the original floating-point noder, which is faster than snap-rounding (but of course has the robustness issues noted above).  Another is to use a special-purpose noder to provide very fast polygonal coverage union.
Union of a polygonal coverage (10x faster with OverlayNG)
  • modular and cleaner codebase allows easier testability, maintenance, enhancement and reuse.  A winged-edge graph model is used for the topology graph. This is simpler and less memory intensive.
  • The rebuild gives an opportunity to make some semantic improvements: 
    • Empty results are returned as empty atomic geometries of appropriate type, rather than awkward-to-handle empty GeometryCollections
    • Linear output is merged node-to-node.  This gives union a more natural and useful  semantic
A benefit of the new codebase is that it is easier to enhance and extend.  For example, it should be straightforward to finally provide a SplitPolygon function for JTS.  Another potential extension is overlay for Polygonal Coverages.

Code that is so widely used needs to be thoroughly tested against real-world workloads.  Initially OverlayNG will be released as a separate API in JTS.  This allows it to be used along with the original overlay.  It can be used as a fallback for cases which fail in the original overlay process.  Once the new code has been proved out in real world use, it is likely to become the standard overlay code path. Also, the code will be ported to GEOS soon, where we're hoping it will provide significant benefits to the many systems that use GEOS.

I'll be posting more articles about aspects of OverlayNG soon. The code is almost ready to release, after some final testing. In the meantime, the pre-release code is available in a Git branch. It would be great to get as much beta-testing as possible before final release, so try it out and log some feedback!

Tuesday, 14 April 2020

Maximum Inscribed Circle and Largest Empty Circle in JTS

There is often a need to find a point which is guaranteed to lie in the interior of a polygon.  Uses include placing cartographic labels, and using the point as a proxy for polygon containment or overlap (such as in polygon overlay).

There are several ways to compute a "centre point" for a polygon.  The simplest is the polygon centroid, which is the Center of Mass of the polygon area.  This has a straightforward O(N) algorithm, but it has the significant downside of not always lying inside the polygon!  (For instance, the centroid of a "U" shape lies outside the shape).  This makes it non-useful as an interior point algorithm.

JTS provides the InteriorPoint algorithm, which is guaranteed to return a point in the interior of a polygon.  This works as follows:
  • Determine a horizontal scan line on which the interior point will be located. To increase the chance of the scan line having non-zero-width intersection with the polygon the scan line Y ordinate is chosen to be near the centre of the polygon's Y extent but distinct from all of vertex Y ordinates.
  • Compute the sections of the scan line which lie in the interior of the polygon.
  • Choose the widest interior section and take its midpoint as the interior point.
This works perfectly for finding proxy points, and usually produces reasonable results as a label point.  

However, there are polygons for which the above algorithm finds points which lie too close to the boundary for a label point.  A better choice would be the point that lies farthest from the edges of the polygon.  The geometric term for this construction is the Maximum Inscribed Circle. The farthest point is the center of the circle. 

Comparison of center points for a not-so-typical polygon

In the geographic domain this is romantically termed the Pole of Inaccessibility.  
Pole Of Inaccessibility in Canada

This point occurs at a node of the medial axis of the polygon, so in theory all that is needed is to compute the medial axis and test the set of node points.  However, medial axis algorithms are notoriously difficult to implement, and can be expensive to compute.  So it's appealing to look for a simple and fast way to compute a good approximation to the Maximum Inner Circle center.  There have been various approaches to this, including a geodetic grid-based approach by Garcia-Castellanos & Lombardo, and one by Martinez using random point distributions.  Recently Mapbox released a clever implementation which uses successive refinement of a grid along with a branch-and-bound technique to reduce the amount of searching needed.  

JTS now has a version of this algorithm, called MaximumInscribedCircle.  It significantly improves performance by using spatial indexing techniques for both polygon interior testing and distance computation.  This makes it very fast to find the MIC even for large, complex polygons.  Performance is key for computing label points, since it is likely to be used for many polygons on a typical map.

Grid refinement to find Maximum Inscribed Circle

An interesting property of the MIC is that its radius is the distance at which the negative buffer of the polygon disappears (becomes empty).  Thus the MIC radius length is a measure of the "narrowness" of a polygon.   This is often useful for purposes of simplification or data cleaning, to remove narrow polygonal artifacts in data.

Sequence of negative buffers containing the Maximum Inscribed Circle center

And, as the infomercials say, that's not all!  If you act today you also get a free implementation of  the Largest Empty Circle !  The Largest Empty Circle is defined for a set of geometric obstacles.  It is the largest circle that can be constructed whose interior does not intersect any obstacle and whose center lies in the convex hull of the obstacles.  The obstacles can be points, lines or polygons (although only the first two are currently implemented in JTS).  Classic use cases for the Largest Empty Circle are in logistics to find a location for a new chain store in a set of store locations; or to find the largest roadless area in environmental planning.

It turns out that the LEC can be computed by essentially the same algorithm as the MIC, with a few small changes.  And of course it also uses spatial indexing to provide excellent performance.
Largest Empty Circles for point and line obstacles

Maximum Inscribed Circle and Largest Empty Circle are now in JTS master, and will be released in the upcoming version 1.17.

Further Improvements

There are some useful enhancements that can be made:
  • For Maximum Inscribed Circle, allow a second polygonal constraint.  This supports finding a label point within a view window rectangle.
  • For Largest Empty Circle, allow a client-defined boundary polygon.  This allows restricting the circle to lie within a tighter bound than the convex hull
  • For both algorithms, it should be feasible to automatically determine a termination tolerance

Monday, 3 February 2020

Running commands in the JTS TestBuilder

The JTS TestBuilder is a great tool for creating geometry, processing it with JTS spatial functions, and visualizing the results.  It has powerful capabilities for inspecting the fine details of geometry (such as the Reveal Topology mode).  I've often thought it would be handy if there was a similar tool for PostGIS.  Of course QGIS excels at visualizing the results of PostGIS queries.  But it doesn't offer the same simplicity for creating geometry and passing it into PostGIS functions.

This is the motivation behind a recent enhancement to the TestBuilder to allow running external (system) commands that return geometry output.  The output can be in any text format that TestBuilder recognizes (currently WKT, WKB and GeoJSON).   It also provides the ability to encode the A and B TestBuilder input geometries as literal WKT or WKB values in the command.  The net result is the ability to run external geometry functions just as if they were functions built into the TestBuilder.


Running PostGIS spatial functions

Combined with the versatile Postgres command-line interface psql, this allows running a SQL statement and loading the output as geometry. Here's an example of running a PostGIS spatial function.  In this case a MultiPoint geometry has been digitized in the TestBuilder, and processed by the ST_VoronoiPolygons function.  The SQL output geometry is displayed as the TestBuilder result.

The command run is:

/Applications/ -qtA -c 
"SELECT ST_VoronoiPolygons('#a#'::geometry);"

Things to note:
  • the full path to psql is needed because the TestBuilder processes the command using a plain sh shell.  (It might be possible to improve this.)
  • The psql options -qtA  suppress messages, table headers, and column alignment, so that only the plain WKB of the result is output
  • The variable #a# has the WKT of the A input geometry substituted when the command is run.  This is converted into a PostGIS geometry via the ::geometry cast.  (#awkb# can be used to supply WKB, if full numeric precision is needed)

Loading data from PostGIS

This also makes it easy to load data from PostGIS to make use of TestBuilder geometry analysis and visualization capabilities.  The query can be any kind of SELECT statement, which makes it easy to control what data is loaded.  For large datasets it can be useful to draw an Area of Interest in the TestBuilder and use that as a spatial filter for the query.  The TestBuilder is able to load multiple textual geometries, so there is not need to collect the query result into a single geometry.

Loading data from the Web

Another use for commands is to load data from the Web, by using curl to retrieve a dataset.  Many web spatial datasets are available in GeoJSON, which loads fine into the TestBuilder. Here's an example of loading a dataset provided by an OGC Features service (pygeoapi):

Command Panel User Interface

The Command panel provides a simple UI to make it easier to work with commands.  Command text  can be pasted and cleared.  A history of commands run is recorded for the TestBuilder session.  Recorded commands in the session can be recalled via the Previous and Next Command buttons.
Buttons are provided to insert substitution variable text.

To help debug incorrect command syntax, error output from commands is displayed.

It can happen that a command executes successfully, but returns output that cannot be parsed.  This is indicated by an error in the Result panel.  A common cause of this is that the command produces logging output as well as the actual geometry text, which interferes with parsing.  To aid in debugging this situation the command window shows the first few hundred characters of command output.  The good news is that many commands offer a "quiet mode" to provide data-only output.

Unparseable psql output due to presence of column headers.  The pqsl -t option fixes this.

If you find an interesting use for the TestBuilder Command capability, post it in the comments!