ST_Subdivide

[Komzpa]

Let ST_Subdivide reselect pivot based on current piece geometry and not the original bbox. Split line for large polygons is also shifted to pass through existing point that lets the subdivision converge faster and chop off larger simple boxes.

This PR passes tests on simple cases.

One of unit tests (about inverted geometries) gets stuck on geos internals:

kom@junocat ~ % sudo gdb -ex "set pagination 0" -ex "thread apply all bt" --batch -p 17978
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
0x00007fe8a2b75d44 in geos::operation::intersection::distance(geos::operation::intersection::Rectangle const&, double, double, double, double) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so

Thread 1 (Thread 0x7fe8b65ff780 (LWP 17978)):
#0  0x00007fe8a2b75d44 in geos::operation::intersection::distance(geos::operation::intersection::Rectangle const&, double, double, double, double) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#1  0x00007fe8a2b76134 in geos::operation::intersection::distance(geos::operation::intersection::Rectangle const&, std::vector<geos::geom::Coordinate, std::allocator<geos::geom::Coordinate> > const&) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#2  0x00007fe8a2b773e5 in geos::operation::intersection::RectangleIntersectionBuilder::reconnectPolygons(geos::operation::intersection::Rectangle const&) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#3  0x00007fe8a2b74c90 in geos::operation::intersection::RectangleIntersection::clip_polygon_to_polygons(geos::geom::Polygon const*, geos::operation::intersection::RectangleIntersectionBuilder&, geos::operation::intersection::Rectangle const&) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#4  0x00007fe8a2b75103 in geos::operation::intersection::RectangleIntersection::clip_geom(geos::geom::Geometry const*, geos::operation::intersection::RectangleIntersectionBuilder&, geos::operation::intersection::Rectangle const&, bool) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#5  0x00007fe8a2b75490 in geos::operation::intersection::RectangleIntersection::clip() () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#6  0x00007fe8a2b75513 in geos::operation::intersection::RectangleIntersection::clip(geos::geom::Geometry const&, geos::operation::intersection::Rectangle const&) () from /usr/lib/x86_64-linux-gnu/libgeos-3.7.0dev.so
#7  0x00007fe8a3c54d30 in GEOSClipByRect_r () from /usr/lib/x86_64-linux-gnu/libgeos_c.so.1
#8  0x00007fe8a3eee2e3 in lwgeom_clip_by_rect (geom1=geom1@entry=0x55bb17b68310, x0=3.0481343214686657e-14, y0=-20000000, x1=x1@entry=20000000, y1=y1@entry=-1) at lwgeom_geos.c:1103
#9  0x00007fe8a3ebe90b in lwgeom_subdivide_recursive (geom=geom@entry=0x55bb17b68310, maxvertices=maxvertices@entry=256, depth=3, depth@entry=2, col=col@entry=0x55bb17b69598, clip=0x55bb17b68360) at lwgeom.c:2396
#10 0x00007fe8a3ebe95f in lwgeom_subdivide_recursive (geom=geom@entry=0x55bb17b69768, maxvertices=maxvertices@entry=256, depth=2, depth@entry=1, col=col@entry=0x55bb17b69598, clip=<optimized out>) at lwgeom.c:2409
#11 0x00007fe8a3ebed0b in lwgeom_subdivide_recursive (geom=geom@entry=0x55bb17b69438, maxvertices=maxvertices@entry=256, depth=1, depth@entry=0, col=col@entry=0x55bb17b69598, clip=clip@entry=0x7fff8826b480) at lwgeom.c:2400
#12 0x00007fe8a3ec0040 in lwgeom_subdivide (geom=0x55bb17b69438, maxvertices=256) at lwgeom.c:2437
#13 0x00007fe8a3e9eb2c in ST_Subdivide (fcinfo=0x55bb17b40eb8) at lwgeom_dump.c:382
#14 0x000055bb1594362c in ExecMakeFunctionResultSet ()
#15 0x000055bb1595c71f in ?? ()
#16 0x000055bb1595c7e4 in ?? ()
#17 0x000055bb159482dc in ?? ()
#18 0x000055bb1594a5e8 in ?? ()
#19 0x000055bb1593c30d in standard_ExecutorRun ()
#20 0x000055bb15a790e6 in ?? ()
#21 0x000055bb15a7a770 in PortalRun ()
#22 0x000055bb15a76193 in ?? ()
#23 0x000055bb15a781c8 in PostgresMain ()
#24 0x000055bb157b1589 in ?? ()
#25 0x000055bb159ffc04 in PostmasterMain ()
#26 0x000055bb157b2e8f in main ()

[pramsey]

So, the pivot is the real vertex closest to the center. And if there is no vertex, it's the center. Yes?

Why is this better/faster? That part I'm not intuiting.

[Komzpa]

@pramsey correct.

It does not create four extra points in case of simple split, but only three, saving one double precision grid snap on diagonal edges, so the polygon can be recombined back with higher precision.

For large polygons with detailed boundary it starts splitting off bigger central 4-point boxes than previously, not insisting on them being squares. After that a balanced pivot should be found in second, now smaller, part faster.

[codecov]

Codecov Report

Merging #202 into svn-trunk will decrease coverage by <.01%.
The diff coverage is 90.47%.

@@              Coverage Diff              @@
##           svn-trunk     #202      +/-   ##
=============================================
- Coverage      79.33%   79.33%   -0.01%     
=============================================
  Files            202      201       -1     
  Lines          63061    63004      -57     
=============================================
- Hits           50032    49982      -50     
+ Misses         13029    13022       -7

Impacted Files	Coverage Δ
liblwgeom/lwutil.c	68.47% <100%> (-0.34%)	⬇️
liblwgeom/lwgeom.c	80.84% <90.16%> (+0.32%)	⬆️
liblwgeom/cunit/cu_lwstroke.c	95.34% <0%> (-1.71%)	⬇️
liblwgeom/lwstroke.c	85.14% <0%> (-0.23%)	⬇️
raster/rt_pg/rtpg_wkb.c
liblwgeom/cunit/cu_ptarray.c	99.79% <0%> (+0.04%)	⬆️
raster/rt_pg/rtpg_band_properties.c	76.34% <0%> (+0.11%)	⬆️
raster/rt_pg/rtpg_inout.c	71.84% <0%> (+4.72%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1d78934...cc5676c. Read the comment docs.

[Komzpa]

So, to get it fixed without bumping geos minimal version we have to fall back on just simple intersection for clipping by box. As a bonus lower versions of geos now support ST_Subdivide too.

[pramsey]

Using full-on intersection doesn't exact a performance penalty?

[mloskot]

Why this deliberately drifts away from the existing coding style?

[dbaston]

We should fail on invalid inputs, not silently run a magic function that may or may not produce a desired result.

[Komzpa]

@pramsey full intersection still finishes in sane time. Either that, or guaranteed endless loop in geos clipbypox2d.

[Komzpa]

On my complex polygon splitting changes from 274466 pieces into 115493 just by reselecting pivot.
Less simple boxes are generated. Number of points in geometry is 206341, so a number of polygons originally is larger than number of input points :)

[Komzpa]

Also, on this geometry some boxes are still lost with ST_Subdivide(geom) by current geos clip by box.