[SCM] PostGIS branch master updated. 3.4.0rc1-792-gd8bbc4f97
git at osgeo.org
git at osgeo.org
Wed Nov 22 13:13:22 PST 2023
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via d8bbc4f9758136ced17ad495463094e093537319 (commit)
from 4e095177d8515b17897bcd70f33ad870746303d4 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit d8bbc4f9758136ced17ad495463094e093537319
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Nov 22 13:12:42 2023 -0800
Improve ST_GeneratePoints performance
for newer GEOS, and in general by
not testing grid squares that are outside of the shape.
Fix memory allocation issue for very narrow inputs.
References #5571
diff --git a/NEWS b/NEWS
index 9d4d58615..6eca201a7 100644
--- a/NEWS
+++ b/NEWS
@@ -20,6 +20,8 @@ To take advantage of all SFCGAL featurs, SFCGAL 1.5.0+ is needed.
- DocBook5 XSL is now required to build html documentation
(Sandro Santilli)
- #5602, Drop support for GEOS 3.6 and 3.7 (Regina Obe)
+ - #5571, Improve ST_GeneratePoints performance, but old
+ seeded pseudo random points will need to be regenerated.
* New Features *
diff --git a/doc/reference_processing.xml b/doc/reference_processing.xml
index 4254b3aae..796f08413 100644
--- a/doc/reference_processing.xml
+++ b/doc/reference_processing.xml
@@ -1277,20 +1277,9 @@ result
<funcsynopsis>
<funcprototype>
<funcdef>geometry <function>ST_GeneratePoints</function></funcdef>
- <paramdef>
- <parameter>g</parameter>
- <type>geometry</type>
- </paramdef>
- <paramdef>
- <parameter>npoints</parameter>
- <type>integer</type>
- </paramdef>
- </funcprototype>
- <funcprototype>
- <funcdef>geometry <function>ST_GeneratePoints</function></funcdef>
- <paramdef> <type>geometry</type> <parameter>g</parameter> </paramdef>
- <paramdef> <type>integer</type> <parameter>npoints</parameter> </paramdef>
- <paramdef> <type>integer</type> <parameter>seed = 0</parameter> </paramdef>
+ <paramdef><type>geometry</type> <parameter>g</parameter></paramdef>
+ <paramdef><type>integer</type> <parameter>npoints</parameter></paramdef>
+ <paramdef choice="opt"><type>integer</type> <parameter>seed = 0</parameter></paramdef>
</funcprototype>
</funcsynopsis>
diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index 44e2187d0..49e5fdf8a 100644
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -1658,11 +1658,11 @@ static void test_point_density(void)
/* Check if the 1000th point is the expected value.
* Note that if the RNG is not portable, this test may fail. */
pt = (LWPOINT*)mpt->geoms[999];
- // ewkt = lwgeom_to_ewkt((LWGEOM*)pt);
+ // char* ewkt = lwgeom_to_ewkt((LWGEOM*)pt);
// printf("%s\n", ewkt);
// lwfree(ewkt);
- CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_x(pt), 0.801167838758, 1e-11);
- CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_y(pt), 0.345281131175, 1e-11);
+ CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_x(pt), 0.363667838758, 1e-11);
+ CU_ASSERT_DOUBLE_EQUAL(lwpoint_get_y(pt), 0.782781131175, 1e-11);
lwmpoint_free(mpt);
pt = NULL;
diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index d546f42ea..d7aa1cedb 100644
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -1399,10 +1399,31 @@ lwgeom_offsetcurve(const LWGEOM* geom, double size, int quadsegs, int joinStyle,
return result;
}
+static GEOSGeometry*
+lwpoly_to_points_makepoly(double xmin, double ymin, double cell_size)
+{
+ GEOSCoordSequence *cs = GEOSCoordSeq_create(5, 2);
+ GEOSCoordSeq_setXY(cs, 0, xmin, ymin);
+ GEOSCoordSeq_setXY(cs, 1, xmin + cell_size, ymin);
+ GEOSCoordSeq_setXY(cs, 2, xmin + cell_size, ymin + cell_size);
+ GEOSCoordSeq_setXY(cs, 3, xmin, ymin + cell_size);
+ GEOSCoordSeq_setXY(cs, 4, xmin, ymin);
+ return GEOSGeom_createPolygon(GEOSGeom_createLinearRing(cs), NULL, 0);
+}
+
LWMPOINT*
lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
{
- double area, bbox_area, bbox_width, bbox_height;
+
+ typedef struct CellCorner {
+ double x;
+ double y;
+ } CellCorner;
+
+ CellCorner* cells;
+ uint32_t num_cells = 0;
+
+ double area, bbox_area, bbox_width, bbox_height, area_ratio;
GBOX bbox;
const LWGEOM* lwgeom = (LWGEOM*)lwpoly;
uint32_t sample_npoints, sample_sqrt, sample_width, sample_height;
@@ -1413,15 +1434,10 @@ lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
uint32_t npoints_tested = 0;
GEOSGeometry* g;
const GEOSPreparedGeometry* gprep;
- GEOSGeometry* gpt;
- GEOSCoordSequence* gseq;
LWMPOINT* mpt;
int32_t srid = lwgeom_get_srid(lwgeom);
int done = 0;
- int* cells;
- const size_t size = 2 * sizeof(int);
- char tmp[2 * sizeof(int)];
- const size_t stride = 2 * sizeof(int);
+
if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
{
@@ -1447,15 +1463,20 @@ lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
return NULL;
}
- /* Gross up our test set a bit to increase odds of getting coverage in one pass */
- sample_npoints = npoints * bbox_area / area;
+ /* Gross up our test set a bit to increase odds of getting coverage
+ * in one pass. For narrow inputs, it is possible we will get
+ * no hits at all. A fall-back approach might use a random stab-line
+ * to find pairs of edges that bound interior area, and generate a point
+ * between those edges, in the case that this gridding system
+ * does not generate the desired number of points */
+ area_ratio = FP_MIN(bbox_area / area, 10000.0);
+ sample_npoints = npoints * area_ratio;
/* We're going to generate points using a sample grid as described
- * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html to try and get a more uniform
- * "random" set of points. So we have to figure out how to stick a grid into our box */
- sample_sqrt = lround(sqrt(sample_npoints));
- if (sample_sqrt == 0)
- sample_sqrt = 1;
+ * http://lin-ear-th-inking.blogspot.ca/2010/05/more-random-points-in-jts.html
+ * to try and get a more uniform "random" set of points. So we have to figure
+ * out how to stick a grid into our box */
+ sample_sqrt = ceil(sqrt(sample_npoints));
/* Calculate the grids we're going to randomize within */
if (bbox_width > bbox_height)
@@ -1481,8 +1502,38 @@ lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
}
gprep = GEOSPrepare(g);
- /* Get an empty multi-point ready to return */
- mpt = lwmpoint_construct_empty(srid, 0, 0);
+
+ /* Now we fill in an array of cells, only using cells that actually
+ * intersect the geometry of interest, and finally weshuffle that array,
+ * so we can visit the cells in random order to avoid visual ugliness
+ * caused by visiting them sequentially */
+ cells = lwalloc(sizeof(CellCorner) * sample_height * sample_width);
+ for (i = 0; i < sample_width; i++)
+ {
+ for (j = 0; j < sample_height; j++)
+ {
+ int intersects = 0;
+ double xmin = i * sample_cell_size;
+ double ymin = j * sample_cell_size;
+ GEOSGeometry *gcell = lwpoly_to_points_makepoly(xmin, ymin, sample_cell_size);
+ intersects = GEOSPreparedIntersects(gprep, gcell);
+ if (intersects == 2)
+ {
+ GEOSGeom_destroy(gcell);
+ GEOSPreparedGeom_destroy(gprep);
+ GEOSGeom_destroy(g);
+ lwerror("%s: GEOS exception on GEOSPreparedIntersects: %s", __func__, lwgeom_geos_errmsg);
+ return NULL;
+ }
+ if (intersects == 1)
+ {
+ cells[num_cells].x = xmin;
+ cells[num_cells].y = ymin;
+ num_cells++;
+ }
+ GEOSGeom_destroy(gcell);
+ }
+ }
/* Initiate random number generator.
* Repeatable numbers are generated with seed values >= 1.
@@ -1490,63 +1541,57 @@ lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
* Unix time (seconds) and process ID. */
lwrandom_set_seed(seed);
- /* Now we fill in an array of cells, and then shuffle that array, */
- /* so we can visit the cells in random order to avoid visual ugliness */
- /* caused by visiting them sequentially */
- cells = lwalloc(2 * sizeof(int) * sample_height * sample_width);
- for (i = 0; i < sample_width; i++)
- {
- for (j = 0; j < sample_height; j++)
- {
- cells[2 * (i * sample_height + j)] = i;
- cells[2 * (i * sample_height + j) + 1] = j;
- }
- }
-
/* Fisher-Yates shuffle */
- n = sample_height * sample_width;
+ n = num_cells;
if (n > 1)
{
for (i = n - 1; i > 0; i--)
{
- size_t j = (size_t)(lwrandom_uniform() * (i + 1));
-
- memcpy(tmp, (char *)cells + j * stride, size);
- memcpy((char *)cells + j * stride, (char *)cells + i * stride, size);
- memcpy((char *)cells + i * stride, tmp, size);
+ CellCorner tmp;
+ j = (uint32_t)(lwrandom_uniform() * (i + 1));
+ memcpy( &tmp, cells + j, sizeof(CellCorner));
+ memcpy(cells + j, cells + i, sizeof(CellCorner));
+ memcpy(cells + i, &tmp, sizeof(CellCorner));
}
}
+ /* Get an empty multi-point ready to return */
+ mpt = lwmpoint_construct_empty(srid, 0, 0);
+
/* Start testing points */
while (npoints_generated < npoints)
{
iterations++;
- for (i = 0; i < sample_width * sample_height; i++)
+ for (i = 0; i < num_cells; i++)
{
int contains = 0;
- double y = bbox.ymin + cells[2 * i] * sample_cell_size;
- double x = bbox.xmin + cells[2 * i + 1] * sample_cell_size;
+ double y = cells[i].y;
+ double x = cells[i].x;
x += lwrandom_uniform() * sample_cell_size;
y += lwrandom_uniform() * sample_cell_size;
if (x >= bbox.xmax || y >= bbox.ymax) continue;
- gseq = GEOSCoordSeq_create(1, 2);
- GEOSCoordSeq_setXY(gseq, 0, x, y);
-
- gpt = GEOSGeom_createPoint(gseq);
-
- contains = GEOSPreparedIntersects(gprep, gpt);
-
- GEOSGeom_destroy(gpt);
+#if POSTGIS_GEOS_VERSION >= 31200
+ contains = GEOSPreparedIntersectsXY(gprep, x, y);
+#else
+ {
+ GEOSGeometry *gpt;
+ GEOSCoordSequence *gseq = GEOSCoordSeq_create(1, 2);;
+ GEOSCoordSeq_setXY(gseq, 0, x, y);
+ gpt = GEOSGeom_createPoint(gseq);
+ contains = GEOSPreparedIntersects(gprep, gpt);
+ GEOSGeom_destroy(gpt);
+ }
+#endif
if (contains == 2)
{
GEOSPreparedGeom_destroy(gprep);
GEOSGeom_destroy(g);
- lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
+ lwerror("%s: GEOS exception on GEOSPreparedIntersects: %s", __func__, lwgeom_geos_errmsg);
return NULL;
}
- if (contains)
+ if (contains == 1)
{
npoints_generated++;
mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
diff --git a/regress/core/regress_brin_index_expected b/regress/core/regress_brin_index_expected
index 5b2c99e3c..2a042a670 100644
--- a/regress/core/regress_brin_index_expected
+++ b/regress/core/regress_brin_index_expected
@@ -1,62 +1,62 @@
scan_seq|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_seq|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_seq|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_seq|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_seq|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_idx|Bitmap Heap Scan,Bitmap Index Scan
2d|1
diff --git a/regress/core/regress_index_expected b/regress/core/regress_index_expected
index 8872c19ef..cfc954ecf 100644
--- a/regress/core/regress_index_expected
+++ b/regress/core/regress_index_expected
@@ -1,21 +1,21 @@
scan_idx|Seq Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
scan_seq|Index Scan
-5323|POINT(134.204197 130.674106)
11208|POINT(126.522745 128.356924)
-27373|POINT(129.481991 125.755641)
-33863|POINT(127.143785 131.932614)
+19845|POINT(127.584643 134.083138)
+27373|POINT(125.017705 130.219927)
+33863|POINT(131.608071 127.468328)
45851|POINT(130.986464 132.890625)
&&|1|5+-5:true
-&&|2|907+-60:true
+&&|2|912+-60:true
&&|3|12505+-500:true
&&|4|50000+-600:true
expr &&|1|5+-5:true
-expr &&|2|907+-60:true
+expr &&|2|912+-60:true
expr &&|3|12505+-500:true
expr &&|4|50000+-600:true
_st_sortablehash|0|768602608280535040|768602608280535040
diff --git a/regress/core/tickets.sql b/regress/core/tickets.sql
index 4f246275d..8fcc885e3 100644
--- a/regress/core/tickets.sql
+++ b/regress/core/tickets.sql
@@ -1252,10 +1252,10 @@ SELECT '#4299', ST_Disjoint(ST_GeneratePoints(g, 1000), ST_GeneratePoints(g, 100
FROM (SELECT 'POLYGON((0 0,1 0,1 1,0 1,0 0))'::geometry AS g) AS f;
-- #4304
-SELECT '#4304', ST_Equals(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000, 12345)),
-ST_Disjoint(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000, 54321)),
-ST_Disjoint(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000)),
-ST_Distance(ST_GeometryN(ST_GeneratePoints(g, 1000, 12345), 1000), 'POINT(0.801167838758 0.345281131175)'::geometry) < 1e-11
+SELECT '#4304',
+ ST_Equals(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000, 12345)),
+ ST_Disjoint(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000, 54321)),
+ ST_Disjoint(ST_GeneratePoints(g, 1000, 12345), ST_GeneratePoints(g, 1000))
FROM (SELECT 'POLYGON((0 0,1 0,1 1,0 1,0 0))'::geometry AS g) AS f;
-- #4331
diff --git a/regress/core/tickets_expected b/regress/core/tickets_expected
index 44e4f6e65..145994491 100644
--- a/regress/core/tickets_expected
+++ b/regress/core/tickets_expected
@@ -384,7 +384,7 @@ equals211|t
equals212|t
equals213|t
#4299|t
-#4304|t|t|t|t
+#4304|t|t|t
ERROR: BOX3D_construct: args can not be empty points
ERROR: BOX2D_construct: args can not be empty points
#4176|t
-----------------------------------------------------------------------
Summary of changes:
NEWS | 2 +
doc/reference_processing.xml | 17 +---
liblwgeom/cunit/cu_algorithm.c | 6 +-
liblwgeom/lwgeom_geos.c | 139 ++++++++++++++++++++-----------
regress/core/regress_brin_index_expected | 60 ++++++-------
regress/core/regress_index_expected | 16 ++--
regress/core/tickets.sql | 8 +-
regress/core/tickets_expected | 2 +-
8 files changed, 143 insertions(+), 107 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list