[PostGIS] #6135: ST_GeneratePoints can return a different number of points than requested for MultiPolygons

[PostGIS]

#6135: ST_GeneratePoints can return a different number of points than requested
for MultiPolygons
---------------------+---------------------------
 Reporter:  Henni    |      Owner:  pramsey
     Type:  defect   |     Status:  new
 Priority:  medium   |  Milestone:  PostGIS 3.6.5
Component:  postgis  |    Version:  3.6.x
 Keywords:           |
---------------------+---------------------------
 Hi,

 I noticed that `ST_GeneratePoints(multipolygon, n)` does not always return
 `n` points, and for `n = 1` it can return NULL. I may be misreading the
 intended behaviour, but the docs suggest exactly `n` points are expected.

 **Version:**
 PostGIS 3.5.3 and 3.6.1 (via imresamu/postgis linux/arm64 docker image)

 **Reproduction** (three equal squares):

 {{{
 WITH g AS (
   SELECT 'MULTIPOLYGON(((0 0,1 0,1 1,0 1,0 0)),
                        ((2 0,3 0,3 1,2 1,2 0)),
                        ((4 0,5 0,5 1,4 1,4 0)))'::geometry AS geom
 )
 SELECT n AS requested, ST_NumGeometries(ST_GeneratePoints(geom, n)) AS
 returned
 FROM g, generate_series(1, 6) AS n;
 }}}

 {{{
 | requested | returned |
 | --------- | -------- |
 |         1 |     NULL |
 |         2 |        3 |
 |         3 |        3 |
 |         4 |        3 |
 |         5 |        6 |
 |         6 |        6 |
 }}}

 The seeded variant `ST_GeneratePoints(geom, n, 42)` behaves the same way.
 A single polygon of the same area returns exactly `n`.

 **Possible cause:**
 This is Claude suggested. I don't have the depth of understanding myself
 to explore this:
 ''`lwmpoly_to_points` in `liblwgeom/lwgeom_geos.c` rounds each part's
 share
 independently (`lround(npoints * sub_area / area)`), so the per-part
 counts don't always add up to `npoints`. This looks unchanged on master.''

 **Possible fix:**
 allocate points to parts with largest-remainder rounding (floor each
 share, then give the leftover
 points to the parts with the largest remainders), so the total is always
 `npoints`.
 I'm happy to help test a patch.

 We ran into this with real municipality boundaries that have exclaves (off
 by ±1 for about 3% of the multi-part cases we tried), and are working
 around it by splitting the MultiPolygon and calling ST_GeneratePoints per
 part.

 Thank you for your work and thanks for PostGIS!
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/6135>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.