Simplify a branching line string

Question

I have been working on developing a GeoTools function to create a central label line for polygons.

The idea was to grab the skeleton of the polygon and use that to place my label, however it turns out that the skeleton of real world polygons is too "spiky" to work as a good label line (labels won't go across the joints).

So I'm looking for a general algorithm that can remove the short dangly edges of this MultiLineString (I also have all the edges in a list if that is easier). Ideally, I'd like Java code based on JTS but I'm happy to work the code up from an algorithm or pseudo code if that is all you have.

Answer 1

A possible algorithm would be to extract the longest path of the skeleton, as a graph structure. This could be obtained by computing all shortest pathes between pairs of extreme leaf nodes (the ones with degree = 1) and keep the longest. (it assumes the skeleton has no cycle, which can happen when the input polygon has holes).

This longest path has great chances to be the line you are looking for label placement, after some smoothing maybe (with that?).

Another possibility could be to use an algorithm relying on the "depth" of a node: Each node would be qualified with its network distance to the nearest leaf. The ridge line of these values could be a good candidate line.

It however does not handle the case when the main line of the skeleton has some more unexpected changes of direction. For that, an algorithm relying on the "stroke" concept presented by:

Thomson, R., 2006. The' stroke' concept in geographic network generalization and analysis. In: Riedl, A., Kainz, W., Elmes, G. (Eds.), Progress in Spatial Data Handling. Springer Berlin Heidelberg, pp. 681-697. URL http://dx.doi.org/10.1007/3-540-35589-8_43

could be interesting to test. And there is luckily an implementation of that in Java/JTS here !

Answer 2

As @Julien suggests (and I saw mentioned on the MapServer PR mentioned in comments by user30184) the best solution I've found is to use Dijkdtra's shortest path algorithm to find the longest shortest path between two nodes of degree 1 in the graph. As far as I can tell this was first proposed by Noah Veltman in his JavaScript code for Centreline placement.

This turns out to be releatively easy to do in GeoTools using the gt-graph module:

  private static Geometry reduceToCentreLine(Geometry geom) {
    LineStringGraphGenerator gen = new LineStringGraphGenerator();
    for (int i = 0; i < geom.getNumGeometries(); i++) {
      gen.add(geom.getGeometryN(i));
    }
    Graph graph = gen.getGraph();
    EdgeWeighter weighter = e -> {
      Geometry g = (Geometry) e.getObject();
      return g.getLength();
    };
    double bestLen = Double.NEGATIVE_INFINITY;
    Path bestPath = null;
    for (Node source : graph.getNodesOfDegree(1)) {
      // calculate the cost(distance) of each graph node to the node closest to
      // the origin
      // System.out.println("starting at " + source.getObject());
      DijkstraShortestPathFinder dspf = new DijkstraShortestPathFinder(graph, source, weighter);
      dspf.calculate();
      for (Node dest : graph.getNodesOfDegree(1)) {
        // System.out.println("\troute to " + dest.getObject());
        if (dest.equals(source)) {
          continue;
        }
        // get length
        double len = 0.0;
        Path path = dspf.getPath(dest);
        if (path == null) {// no path to dest
          continue;
        }
        for (Edge e : path.getEdges()) {
          Geometry g = (Geometry) e.getObject();
          len += g.getLength();
        }
        if (len > bestLen) {
          bestPath = path;
          bestLen = len;
        }
      }
    }
    ArrayList<LineString> edges = new ArrayList<>();
    for (Edge e : bestPath.getEdges()) {
      Geometry g = (Geometry) e.getObject();
      edges.add((LineString) g);
    }
    return GF.createMultiLineString(GeometryFactory.toLineStringArray(edges));
  }

There is still some work to do with the densification and simplification steps and some speed up but I'm quite happy with the results so far:

Answer 3

So, I am publishing my new experimental improved function ST_CentraLAxisLongestLineFromVoronoiDiagrams, to try to make you a little happier :-), and for JAVA programmers to try to implement this approach based on the Voronoi Diagrams function.

Create a function:

CREATE OR REPLACE FUNCTION ST_CentraLAxisLongestLineFromVoronoiDiagrams(
geom GEOMETRY,
n integer)
RETURNS GEOMETRY AS
$BODY$
WITH
    tbla AS (SELECT (ST_Dump(geom)).geom),
    tblb AS (SELECT ST_ExteriorRing(geom) geom1, ST_Buffer(ST_ExteriorRing(geom), 0.000001) geom2 FROM tbla),
    tblc AS (SELECT line1, line2, line3 FROM (SELECT ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),1)),
(ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),2))) AS line1, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),3)), 
(ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),4))) AS line2, ST_MakeLine((ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),1)), 
(ST_PointN(ST_Boundary(ST_OrientedEnvelope(geom1)),4))) AS line3 FROM tblb) AS moo),
    tbld AS (SELECT (ST_Dump(ST_Difference(a.geom2, ST_Union(b.geom)))).geom AS geom FROM tblb a CROSS JOIN LATERAL
             (SELECT ST_Buffer(line1, 0.00001) geom FROM tblc UNION SELECT ST_Buffer(line2, 0.00001) geom FROM tblc) b GROUP BY a.geom2),
    tble AS (SELECT (a.geom) geom FROM tbld a JOIN tblc b ON ST_Intersects(a.geom, b.line3)),
    tblf AS (SELECT generate_series (0, n) as steps),
    tblg AS (SELECT steps AS stp, ST_LineInterpolatePoint(geom1, steps/(SELECT count(steps)::float-1 FROM tblf)) geom FROM tblb, tblf GROUP BY tblf.steps, geom),
    tblh AS (SELECT ((ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom) geom FROM tblg),
    tbli AS (SELECT ST_ExteriorRing(ST_Union(a.geom)) geom FROM tblh a JOIN tble b ON ST_Intersects(a.geom, b.geom))
            SELECT ST_Intersection(a.geom, b.geom) geom FROM tbli a JOIN tbla b ON ST_Intersects(a.geom, b.geom);
$BODY$
LANGUAGE SQL;

Run the function on polygons and check the result.

SELECT ST_CentraLAxisLongestLineFromVoronoiDiagrams(geom, 789) geom FROM <polygon_table>

To improve the result, you can apply ST_Simplify(), as increasing the number of points also improves the result, but negatively affects the performance...

For example:

SELECT ST_Simplify(geom, 0.01) geom FROM (SELECT ST_CentraLAxisLongestLineFromVoronoiDiagrams(geom, 789) geom FROM <polygon_table>) foo

You may need to adjust the operation of the geo-tool, with ST_Buffer().

I've been playing around with the data in WGS84.

The main thing is all...

Original spatial solutions :-)