lwgeom_mindistance2d slows down topology loading

[strk]

Profiling an slow addition of a single point I found the CPU busy in computing the minimum distance between the newly added edge with 43520 vertices and each of 2862 other edges intersecting it. Such operation takes a very long time, even multiple seconds for each single distance computation.

Overall, the addition of a that _single_ line takes over 13 minutes.

The same long timings can be seen using ST_Distance itself. Examples:

LOG:  duration: 1043.816 ms  statement: select a.id aid, b.id bid, st_npoints(a.wkb_geometry), st_npoints(b.wkb_geometry), st_distance(a.wkb_geometry, b.wkb_geometry) from contour_onet a, contour_onet b where a.id = 14408 and b.id = 9928;
  aid  | bid  | st_npoints | st_npoints |   st_distance    
-------+------+------------+------------+------------------
 14408 | 9928 |      43520 |        499 | 10568.8576870335
(1 row)

LOG:  duration: 4937.063 ms  statement: select a.id aid, b.id bid, st_npoints(a.wkb_geometry), st_npoints(b.wkb_geometry), st_distance(a.wkb_geometry, b.wkb_geometry) from contour_onet a, contour_onet b where a.id = 14408 and st_npoints(b.wkb_geometry) = 2247;
  aid  | bid  | st_npoints | st_npoints |   st_distance
-------+------+------------+------------+------------------
 14408 | 9031 |      43520 |       2247 | 53917.1247457271
(1 row)

This ticket is to find speed improvement opportunities.

[strk]

Another example:

LOG:  duration: 6416.539 ms  statement: select st_summary(a), st_summary(b), st_distance(a,b) from slowdist;
            st_summary            |           st_summary            |   st_distance   
----------------------------------+---------------------------------+-----------------
 LineString[BS] with 43520 points | LineString[BS] with 3153 points | 254918.45088224
(1 row)

[strk]

In 14977:

Make adding a line to topology interruptible

… as it otherwise gets frustrating when slow (see #3587)

[nicklas]

Yes, this is good.

I think Paul have the solution, we just have to encourage him a little.

Some background: My first contribution to PostGIS was about this, added in 1.5.

But the result of that only helps when bounding boxes is not overlapping. What I did is described in: ​https://trac.osgeo.org/postgis/wiki/NewDistCalcGeom2Geom

This gives a good boost, when bbox isn't overlapping.

I added 3 geometries as example in attachement.

I copied and moved one of them so the bounding boxes overlap.

So those 2 queries is using the same geometries, but the first query with overlapping bounding boxes:

select st_npoints(a.geom), st_npoints(b.geom), st_distance(a.geom, b.geom) from 
(select geom from t where gid = 1) a,
(select geom from t where gid = 2) b

st_npoints | st_npoints |   st_distance   
------------+------------+-----------------
      55773 |      31248 | 1791.9837591391
(1 row)


Time: 73239.353 ms

select st_npoints(a.geom), st_npoints(b.geom), st_distance(a.geom, b.geom) from 
(select geom from t where gid = 1) a,
(select geom from t where gid = 3) b

 st_npoints | st_npoints |   st_distance    
------------+------------+------------------
      55773 |      31248 | 194387.282408659
(1 row)

Time: 709.125 ms

So this is good, but useless in many situations. But Paul have some tree-structure that will do the thing.

Do it Paul, just do it

[strk]

Is there any way to access the older algorithm, for test ? I think the tree structure you are talking about is the one in postgis/ ? Or is that one for the new algorithm you are talking about ? Whatever it is, having access to that tree in liblwgeom would greatly help…

[nicklas]

That tree thing, I think Paul have implemented in geography, but I am no sure if it is the same thing he talked about to implement in geometry 5 years ago.

About the old algorithm, I don't get what you mean. In 1.4 and before all distance calculations was done eith a "brute force" method. Just comparing all combinations of vertex-vertex and vertex-edge. In 1.5 I implemented a more efficient method for cases with linestring and polygons involved, and bboxes that don't intersect. For cases with overlapping bboxes it falls back to the same old "brute force" method, so it is still there. Not exectly the same as in 1.4, since I did some smaller changes. Since 1.5 I think the distance code for geometry is more or less is untouched except that someone (Paul I think) have added support for circular geometries.

[strk]

In my case I'm only calling lwgeom_mindistance2d against linestring pairs with intersecting bounding boxes. So you're saying I do am seeing the brute force in action ? No other overhead other than a re-check for bounding box overlap ?

Paul: is the rtree cache used for ST_Distance(geometry, geometry) ? Any plan to move the rtree interface down to liblwgeom for use from the topology code ? (it is common to test the same long line against hundreds of other lines, so a local temporary cache could be very useful)

[pramsey]

There is no rtree cache used for distance in geometry currently. The start of the code is in liblwgeom/lwtree.c but it has never been hooked up and is probably rotten by now. I got pulled into circular geometry distance by the requirements of the tree, and never got back to it. The tree code in geography turned up enough cases where brute_distance() != tree_distance() that I have viewed hooking up the tree code as a huge task requiring piles of testing and and QA/QC. It would certainly be very useful for repeated distance testing, where the repeated side is complex.

The only place where we have tree code that has not been pulled down into liblwgeom is the tree code for point-in-poly intersects true/false testing.

[nicklas]

Yes, you are meeting the brute force wall at once the bboxes is overlapping.

About the overhead, things can probably be optimized. But my guess is that will not give more than maybe 10 % increase of speed. If I recall right I had the brute force part run approx 10 % faster when I rewrote it for handling distance and shortest line the same way. The purpose was to get exactly the same distance as the length of the shortest line. As a bonus there was a small increase in speed.

[strk]

I'm still surprised by the excess time. Note that using ST_Intersects(ST_Buffer(a,10), b) takes about 200ms. Given the topology use for mindistance2d is to find out if the geometry is within a given amount of units from each other, I guess the check could be replaced with those two calls to GEOS.

But the extreme difference in time seems to suggest there's something fundamentally wrong with the distance algorithm. Is there any way to access the older algorithm, fully brute force, from SQL ?

[strk]

See the numbers:

LOG:  duration: 6360.923 ms  statement: select st_dwithin(a,b,10) from slowdist;
 st_dwithin 
------------
 f
LOG:  duration: 210.333 ms  statement: select st_intersects(st_buffer(a,10),b) from slowdist;
 st_intersects 
---------------
 f
LOG:  duration: 59.426 ms  statement: select st_intersects(a,st_buffer(b,10)) from slowdist;
 st_intersects 
---------------
 f

[nicklas]

Yes, I agree it looks very strange. If you want to test the version before my contribution you should fire up a 1.4.

It should be fairly easy to force the brute force method in todays code too. I can take a look where to search. But that will not rewind to what it was before I added things.

I will try to take a look too if I can see what is happening.

Can the prepared geometry help the intersects case, or is it only when quering many geometries that help?

[strk]

The numbers you see are w/out prepared geometries, I think, because a single ST_Intersects call is issued within the statement.

[strk]

PostGIS 1.4 is actually worst:

pg14=# select st_intersects(a,st_buffer(b,10)) from slowdist;
 st_intersects 
---------------
 f
(1 row)

Time: 69.658 ms
pg14=# select st_intersects(st_buffer(a,10),b) from slowdist;
 st_intersects 
---------------
 f
(1 row)

Time: 423.872 ms
pg14=# select st_dwithin(a,b,10) from slowdist;
 st_dwithin 
------------
 f
(1 row)

Time: 13152.853 ms

[nicklas]

A patch that comments out the code that performs the bbox test and sends to faster distance calc

[nicklas]

I have attached a patch that comments out the bbox test and sending to faster distance calc.

I see no difference in timing with that.

I also think this is very strange.

I can see 3 possible reasons

1) Something very smart is happening in GEOS in st_intersects code

2) We have some very big overhead in our helping functions that for instance do very many memory allocations or reallocations, that I don't realize when I look at the code.

3) There is some ugly brain fart in the distance code itself. But since 1.4 shows the same pattern, I guess that is something we have had for a very long time and that I failed to see when I did my changes.

An interesting thing is that ST_3DIntersects behaves in about the same way as the distance functions. That is not strange since it is a wrapper around similar code as the 2D distance functions.

[pramsey]

My bet is on (1), as I recall a conversation w/ Martin in which is said his distance implementation used a sweepline implementation. He may be tossing all the edges into a bin and discarding many of them far more efficiently than we are. Most of your enhancements were around making multipart objects more efficient, but I think strk's examples are basically one huge part, right? It's possible my tree implementation would be as good as Martin's implementation, but it might not necessarily, as a sweepline could get a more consistent result than a depth-first search, which I'm pretty sure is what I'm getting out of my recursing on the trees.

[strk]

A dump of the "slowdist" subject table

[strk]

I've attached the "slodist" table used for my times above, so you can reproduce. I also bet on GEOS being smart, but note this is a relate operation, not a distance operation. You get the same fast times with ST_Relate(ST_Buffer(..))

[strk]

I confirm objects are single-part:

• LineString[BS] with 43520 points
• LineString[BS] with 3153 points

[nicklas]

small correction @pramsey, the faster way of doing it is all about single geometries. But I couldn't fins a way of using it when geometries is too close.

Interesting also that when I tried it on your geoemtries @strk, it really showed that it doesn't do a goog work when such a long geometry as a is close to the smaller one.

I tried by translating the geometry away from the bounding box to make the faster algorithm kick in. When I move it just outsinde the bounding box, the timing was even worse than the brute force version. But when I move it further away it start to kick ass.

without moving anything on my box:

slowdist=# select st_distance(a, b) from slowdist;
   st_distance   
-----------------
 254918.45088224
(1 row)

Time: 7134.470 ms

Just outside the bbox (worst case inside the "circle" of the long line):

slowdist=# select st_distance(a, st_translate(b,-290000,0)) from slowdist;
   st_distance    
------------------
 432666.185703918
(1 row)


Time: 9459.979 ms

But if I move the geometries more apart so geometry a is more in one direction from b (happens easier with geoemtries that is not so long as a) then we see the difference:

slowdist=# select st_distance(a, st_translate(b,-29000000,0)) from slowdist;
   st_distance    
------------------
 28728819.1262113
(1 row)

Time: 23.064 ms

I haven't seen before cases when the "faster" way is slower than brute force, but I guess this is a little extreme case.

If we insted move geometry b just above the bounding box of a we get nicer results at once:

slowdist=# select st_distance(a, st_translate(b,0,1950000)) from slowdist;
   st_distance    
------------------
 79242.5337417284
(1 row)

Time: 250.828 ms

This doesn't help the case.

But if @pramsey means that geos have also distance algorithms with the performance of ST_Intersects, I guess we should switch to them at once.

If geos (jts) is that good also at distance calculations I think that is something Martin has implemented in the last years. I remember Martin reviewing my work and mentioned that he had thought in the same direction but not tested it. He also confirmed that it will not work on too close geometeries.

[nicklas]

In strks example we get 43520 * 3153 combinations to test. For all those combinations we have to test both point to point and point to segment.

So I guess it is possible that our computers isn't faster than this. If using 7 seconds for the distance caclulation it means about 19000 combinations per millisecond. 19000 point to point and 19000 point to segment (well I guess it is two less point to segment combination).

But there might be things to do about it also with brute force. I mean if we can decrease with 1 microsecond per combination it will give a huge effect.

When doing it the brute force way, I think it is 2 places to look for optimization.

1) Both for point/point and point/segment tests we end up with a sqrt calculation. When I worked on those things long time ago I tested to not do that sqrt calc and only do it on the resulting closest points. But I ran into many strange memory problems. I didn't investigate to the ens, but I guess it results in handling very big numbers which double precision can't handle. I mean 10000 km results in handling the number 10 raised by 14. But if that is solvable (by me just doing something wrong back then) I guess that could make a noticable difference to avoid 38000 sqrt calculations per millisecond.

2) at line 2075 in measures.c is a calculation with a lot of addition, substraction, multiplication and division. That line is executed for each point/segment test. There is more below in the same function, but that is only executed when the projected point on the segment is actually on the segment. More common is that the projected point is outside the segment, and the function goes directly to compare the endpoint of the segment with the point.

[strk]

There's another tree implementation in lwgeom_geos_cluster.c, simply wrapping the GEOS STRTree class. In case we want to implement something like that internally. Or GEOS C-API does have a GEOSDistance function if we want to try that.

About avoiding the sqrt, the geometries could be shifted to the origin to lower the numbers

[strk]

FYI: the GEOSDistance implementation is slower than the PostGIS one, for this case: 30 seconds vs. 15 seconds (running on battery).

For sure it's not the same code path issued by ST_Intersects…

[strk]

The offending lines: A is the long one, B the small one

[strk]

I confirm the MonotoneChain / Sweepline based algorithm is used for constructing the lines topologies (for ST_Intersects / ST_Relate).

If we take the STRTree route, GEOS-3.6.0 (unreleased) will have a GEOSSTRtree_nearest function, so we could build a tree containing the envelopes of each segment (or monotonechain) of the longest geometry and then query it for closest segment (or chain) on each of the shorter geometry segments.

[pramsey]

No, you don't want to build an strtree of the edges, as that will take quite a bit of computational time itself. The trick to trees is to use the natural autocorrelation of the edges to group the leaves into internal nodes. That way no sorting or individual insertion is required, you just roll up the ordered set of edges into a tree in O(N) time, then you recurse back down the tree to find the answer in O(log(N)) time, so in theory the whole thing is never more than O(N) + O(log(N)) = O(N). Now, why Intersects() is still quite fast remains a mystery.

[strk]

I've just tried creating the STRtree and it does not seem to affect the timings much. No querying, just creation and insert. The time is still ~15 seconds.

[strk]

querying the STRtree for nearest segment gets to the distance result in ~20ms. Looks promosing. If only the result was similar to the original one …

LOG:  duration: 19.001 ms  statement: select st_summary(a), st_summary(b), st_distance(a,b) from slowdist;
            st_summary            |           st_summary            |   st_distance
----------------------------------+---------------------------------+-----------------
 LineString[BS] with 43520 points | LineString[BS] with 3153 points | 827058.14279572
(1 row)

[strk]

Ok, a correct implementation takes ~2.3 seconds. It's not 20ms but still much faster than ~13 seconds…

LOG:  duration: 2320.946 ms  statement: select st_distance(b,a) from slowdist;
   st_distance   
-----------------
 254918.45088224
(1 row)
LOG:  duration: 2330.856 ms  statement: select st_dwithin(b,a,10) from slowdist;
 st_dwithin 
------------
 f
(1 row)

It just needs to be cleaned up to go.

Pull request follows:

The following changes since commit f4d8771b218fcf2fc66db3f1ffeecda3b4e8a140:

  Remove trailing spaces (2016-07-04 23:50:26 +0200)

are available in the git repository at:

  https://git.osgeo.org/gogs/postgis/postgis.git svn-trunk-strtree-mindist

for you to fetch changes up to c1e163fa3883b75094623cdd66335f9d4d973727:

  Use STRtree based mindistance computation (2016-07-05 01:12:20 +0200)

----------------------------------------------------------------
Sandro Santilli (1):
      Use STRtree based mindistance computation

 liblwgeom/measures.c | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 61 insertions(+), 1 deletion(-)

Tests are welcome!

Regina, can we find a way to point Jenkins at Gogs pull requests, maybe ?

[strk]

Daniel, given I'm getting a segment which doesn't seem to be the closest one from GEOSSTRtree_nearest, maybe you want to take a look at the PR too ? There could be a bug in GEOS.

What I'm doing:

1. Creating an STRTree with node capacity 10
2. Populating the tree with 43519 segments
3. Finding the closest segment to each of 3152 segments of the other geometry
4. Computing the closest points within the segment pair

The problem is that the brute force approach finds a pair of segments which are closer than the ones found by GEOSSTRtree_nearest

[strk]

Daniel: forget it, it was an error in my code, now works correctly, ad it is still 15 times faster than before:

LOG:  duration: 954.569 ms  statement: select ST_AsText(st_shortestline(a,b)) from slowdist;
                              st_astext                              
---------------------------------------------------------------------
 LINESTRING(5161330.6 -4881460.77503163,5046321.13099501 -4653960.8)
(1 row)

Not as good as ST_Intersection, but still better than what we have now. I guess some experimenting could be done with different nodes capacity (10 sounds a bit low)

[strk]

An STRtree node capacity of 256 resulted in raising the timings to ~5 seconds. A node capacity of 32 made timings =~ 2 seconds. A node capacity of 4 made timings =~ 686 ms. A node capacity of 2 made timings =~ 532 ms. A node capacity of 1 resulted in a segfault.

I'm not sure what to make of these numbers.

[dbaston]

If you're willing to require GEOS 3.6, we should probably port IndexedFacetDistance over from JTS. From a quick glance, it looks like all of the non-trivial code was already captured with the MinimumClearance port. This would avoid needing to build a GEOS LineString for every line segment.

[strk]

Good point about the LineString-per-segment need removal. I've to say IndexedFacetDistance also already has a good API for use from the topology code, so good idea !

[strk]

Another great idea (for topology): DROP the distance computation. Given that distance is only used to reduce the number of ST_Snap operations, and that ST_Snap would be a no-op if no vertex is closer than the given tolerated distance… we might as well just avoid to perform the pre-flight check.

The testsuite is fine with that, which is a good sign!

[dbaston]

FWIW, I ported IndexedFacetDistance to GEOS and exposed it in the CAPI as GEOSDistanceIndexed. Using the two geometries from the table strk posted, the calculation timings are as follows (includes overhead of reading text file and parsing WKT)

GEOSDistance: 41.7 s

GEOSDistanceIndexed: 0.084s

This might not be the right fit for the situation described in this ticket, but it seems like it could be useful to PostGIS overall.

[strk]

I've revived my old pull request on Gitea: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/38

We have CI running for PRs on Gitea now, so no need for Jenkins to do that.

It's to be checked against GEOSDistanceIndexed as reported by dbaston

[strk]

For the record: GEOS-3.9 will have Distance support for PreparedGeometry objects. This would allow building the index only once for the big line. It should be considered as another alternative.

[strk]

The PreparedGeometry approach is in ​https://gitlab.com/postgis/postgis/-/merge_requests/24

[strk]

This is still a problem hit downstream, as shown by ​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/issues/73

[strk]

For the case described in ​https://gitlab.com/nibioopensource/resolve-overlap-and-gap/-/issues/73 we have ST_Intersects(ST_Buffer()) taking 70ms instead of 23 seconds. I don't know if it will be always faster but sounds like something worth testing.

[strk]

WIP PR to use ST_Intersects(ST_Buffer()) is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/179

[strk]

GEOS-3.10.0 introduced GEOSDistanceWithin - performance of using it as as good as the ST_Intersects(ST_Buffer()) one for this case, with no downside of needing to construct the buffer, with the related robustness issue. PR is here: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/180

[strk]

See #2614 for making ST_DWithin hook to GEOS

[strk]

I've attached test10input_singlepoly.sql.gz​ containing the linestrings forming a big polygon. Adding each line with TopoGeo_addLinestring, ordered by "path" (exterior ring first, interior ring later) have the following times on my laptop:

==> load2-bigpoly.log.dist2d_fast_no_check_overlap <==
Time: 1224127.033 ms (20:24.127)

==> load2-bigpoly.log.main <==
Time: 1045277.101 ms (17:25.277)

==> load2-bigpoly.log.topo-addline-geosdistanceindexed <==
Time: 262539.817 ms (04:22.540)

==> load2-bigpoly.log.topo-addline-intersects-buffer <==
Time: 275923.548 ms (04:35.924)

==> load2-bigpoly.log.topo-adline-dwithin <==
Time: 260558.026 ms (04:20.558)

To recap: 17m with master branch, 20m when forcing tree based mindist2d computation, ~4m when using either of the GEOS versions: intersects(buffer()), distanceIndexed, distanceWithin

[strk]

The DistanceIndexed function of GEOS is available since version 3.7 while we require 3.8 in master branch so it would come at no cost to use that

[strk]

It's useful to note that subdividing the lines before adding them into the topology has a great benefit in term of load time. Using ST_Subdivide to split them to have at most 256 vertices makes the whole data load in under 30 seconds, with no change in the PostGIS code (so it would compare to the 17 minutes run).

[robe]

Wow that is an impressive improvement.

[strk]

Yes the speed improvement is impressive but there are two downsides:

1. Current ST_Subdivide introduces new vertices (preventing which was the topic of #4683)
2. The resulting topology has more nodes and edges than strictly required (not necessarely a bad thing)

[Sandro Santilli <strk@…>]

In 865bfe8/git:

Add a topology performance test based on concentric linear rings

References #3587

[strk]

These are the timings of various approaches with the new performance testcase "concentric circles" (adding lines forming concentric circles with 128 segments per quadrant each):

TopoGeo_addLinestring_concentricCircles.log.main
Time: 31703.711 ms (00:31.704)
Time: 31486.980 ms (00:31.487)
Time: 31618.496 ms (00:31.618)
TopoGeo_addLinestring_concentricCircles.log.topo-addline-dwithin
Time: 6832.686 ms (00:06.833)
Time: 6729.710 ms (00:06.730)
Time: 6724.133 ms (00:06.724)
TopoGeo_addLinestring_concentricCircles.log.topo-addline-geosdistanceindexed
Time: 4084.981 ms (00:04.085)
Time: 4071.703 ms (00:04.072)
Time: 4082.751 ms (00:04.083)

It looks like GEOSDistanceIndexed is the best choice as it doesn't change semantic of the code, is the fastest in this case and is available in all supported versions of GEOS.

The timings of the simple case of full coverage are not affected:

TopoGeo_addLinestring.log.main
Time: 44484.710 ms (00:44.485)
Time: 44450.532 ms (00:44.451)
TopoGeo_addLinestring.log.topo-addline-dwithin
Time: 44568.118 ms (00:44.568)
Time: 44238.656 ms (00:44.239)
TopoGeo_addLinestring.log.topo-addline-geosdistanceindexed
Time: 44405.012 ms (00:44.405)
Time: 44752.527 ms (00:44.753)

So the winner is: ​https://git.osgeo.org/gitea/postgis/postgis/pulls/182

[Sandro Santilli <strk@…>]

In 6a6d7da/git:

Use GEOS-3.7+ GEOSDistanceIndexed to find existing edges

Closes #3587 in master branch (3.5.0dev)

[strk]

For the record: further improvement in ​https://git.osgeo.org/gitea/postgis/postgis/pulls/192