- #4272, Improve notice message when unable to compute stats (Raúl Marín)
- #4314, ST_ClipByBox2D: Do not throw when the geometry is invalid (Raúl Marín)
- #4313, #4307, PostgreSQL 12 compatibility (Laurenz Albe, Raúl Marín)
+ - #4299, #4304, ST_GeneratePoints is now VOLATILE. IMMUTABLE version with
+ seed parameter added. (Mike Taves)
PostGIS 2.5.0
Style2;POLYGON((142.928932188135 157.071067811865,144.444297669804 158.314696123025,146.173165676349 159.238795325113,148.049096779839 159.807852804032,150 160,151.950903220161 159.807852804032,153.826834323651 159.238795325113,155.555702330196 158.314696123025,157.071067811865 157.071067811865,158.314696123025 155.555702330196,159.238795325113 153.826834323651,159.807852804032 151.950903220161,160 150,160 50,159.807852804032 48.0490967798387,159.238795325113 46.1731656763491,158.314696123025 44.444297669804,157.071067811865 42.9289321881345,155.555702330196 41.6853038769745,153.826834323651 40.7612046748871,151.950903220161 40.1921471959677,150 40,148.049096779839 40.1921471959677,146.173165676349 40.7612046748871,144.444297669804 41.6853038769745,142.928932188135 42.9289321881345,141.685303876975 44.444297669804,140.761204674887 46.1731656763491,140.192147195968 48.0490967798387,140 50,140 125.857864376269,57.0710678118655 42.9289321881345,55.555702330196 41.6853038769745,53.8268343236509 40.7612046748871,51.9509032201613 40.1921471959677,50 40,48.0490967798387 40.1921471959677,46.1731656763491 40.7612046748871,44.444297669804 41.6853038769745,42.9289321881345 42.9289321881345,41.6853038769746 44.444297669804,40.7612046748871 46.1731656763491,40.1921471959677 48.0490967798387,40 50,40.1921471959677 51.9509032201613,40.7612046748871 53.8268343236509,41.6853038769745 55.555702330196,42.9289321881345 57.0710678118655,142.928932188135 157.071067811865))\r
-Style1-thinline;MULTIPOINT(102.722861415448 109.63774529252,149.121982482376 109.73967711417,143.583483382672 71.2656025879696,43.6127811517685 46.1928159428694,64.0937528611103 64.3971068453017,158.682821131016 128.593401898251,144.217047639393 142.958464308603,143.193456831568 81.4453566087832,137.923520615253 135.276955473495,149.037141026032 48.2552568132572,85.2284310434278 92.0542008728294,107.117526779992 101.481368449965)\r
+Style1-thinline;MULTIPOINT(124.719467098431 131.257475277821,155.60408218128 72.9518237434817,106.69719046413 119.774229927365,140.578248132556 153.822320976712,152.641957278627 87.2817137366839,143.433601275019 137.83621232774,94.5170191740369 89.6648113250309,157.562215520436 47.0045167465619,51.5256023126078 40.4134116206039,158.800069139342 117.796996865768,67.1788241389208 69.3638368956401,129.504207562589 133.236854414052)\r
</paramdef>
<paramdef>
<parameter>npoints</parameter>
- <type>numeric</type>
+ <type>integer</type>
+ </paramdef>
+ </funcprototype>
+ <funcprototype>
+ <funcdef>geometry <function>ST_GeneratePoints</function></funcdef>
+ <paramdef>
+ <parameter>g</parameter>
+ <type>geometry</type>
+ </paramdef>
+ <paramdef>
+ <parameter>npoints</parameter>
+ <type>integer</type>
+ </paramdef>
+ <paramdef>
+ <parameter>seed</parameter>
+ <type>integer</type>
</paramdef>
</funcprototype>
<para>
ST_GeneratePoints generates pseudo-random points until the requested number are
- found within the input area.
+ found within the input area. The optional seed must be greater than zero,
+ and is used to regenerate a deterministic sequence of points.
</para>
<para>Availability: 2.3.0</para>
+ <para>Enhanced: 3.0.0, added seed parameter</para>
</refsection>
<refsection>
<imageobject>
<imagedata fileref="images/st_generatepoints02.png" />
</imageobject>
- <caption><para>Generated 12 Points overlaid on top of original polygon</para></caption>
+ <caption><para>Generated 12 Points overlaid on top of original polygon using a random seed value 1996</para></caption>
</mediaobject>
</informalfigure>
- <programlisting>SELECT ST_GeneratePoints(
- ST_Buffer(
+ <programlisting>SELECT ST_GeneratePoints(geom, 12, 1996)
+FROM (
+ SELECT ST_Buffer(
ST_GeomFromText(
- 'LINESTRING(50 50,150 150,150 50)'
- ), 10, 'endcap=round join=round'), 12);</programlisting>
+ 'LINESTRING(50 50,150 150,150 50)'),
+ 10, 'endcap=round join=round') AS geom
+) AS s;</programlisting>
</para></entry>
</row>
</tbody>
g_util.o \
lwgeodetic.o \
lwgeodetic_tree.o \
+ lwrandom.o \
lwtree.o \
lwout_gml.o \
lwout_kml.o \
lwin_wkt.h \
lwin_wkt_parse.h \
lwout_twkb.h \
+ lwrandom.h \
lwtree.h \
measures3d.h \
measures.h \
{
LWGEOM *geom;
LWMPOINT *mpt;
+ LWMPOINT *mpt2;
+ LWPOINT *pt;
+ LWPOINT *pt2;
+ int eq, i;
// char *ewkt;
/* POLYGON */
- geom = lwgeom_from_wkt("POLYGON((1 0,0 1,1 2,2 1,1 0))", LW_PARSER_CHECK_NONE);
- mpt = lwgeom_to_points(geom, 100);
+ geom = lwgeom_from_wkt("POLYGON((0 0,1 0,1 1,0 1,0 0))", LW_PARSER_CHECK_NONE);
+ mpt = lwgeom_to_points(geom, 100, 0); /* Set a zero seed to base it on Unix time and process ID */
CU_ASSERT_EQUAL(mpt->ngeoms,100);
- // ewkt = lwgeom_to_ewkt((LWGEOM*)mpt);
- // printf("%s\n", ewkt);
- // lwfree(ewkt);
+
+ /* Run a second time with a zero seed to get a different multipoint sequence */
+ mpt2 = lwgeom_to_points(geom, 100, 0);
+ eq = 0;
+ for (i = 0; i < 100; i++)
+ {
+ pt = (LWPOINT*)mpt->geoms[i];
+ pt2 = (LWPOINT*)mpt2->geoms[i];
+ if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
+ eq++;
+ }
+ CU_ASSERT_EQUAL(eq, 0);
lwmpoint_free(mpt);
+ lwmpoint_free(mpt2);
+ pt = NULL;
+ pt2 = NULL;
- mpt = lwgeom_to_points(geom, 1);
- CU_ASSERT_EQUAL(mpt->ngeoms,1);
+ /* Set seed to get a deterministic sequence */
+ mpt = lwgeom_to_points(geom, 1000, 12345);
+
+ /* Check to find a different multipoint sequence with different seed */
+ mpt2 = lwgeom_to_points(geom, 1000, 54321);
+ eq = 0;
+ for (i = 0; i < 1000; i++)
+ {
+ pt = (LWPOINT*)mpt->geoms[i];
+ pt2 = (LWPOINT*)mpt2->geoms[i];
+ if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
+ eq++;
+ }
+ CU_ASSERT_EQUAL(eq, 0);
+ lwmpoint_free(mpt2);
+ pt = NULL;
+ pt2 = NULL;
+
+ /* Check to find an identical multipoint sequence with same seed */
+ mpt2 = lwgeom_to_points(geom, 1000, 12345);
+ eq = 0;
+ for (i = 0; i < 1000; i++)
+ {
+ pt = (LWPOINT*)mpt->geoms[i];
+ pt2 = (LWPOINT*)mpt2->geoms[i];
+ if (lwpoint_get_x(pt) == lwpoint_get_x(pt2) && lwpoint_get_y(pt) == lwpoint_get_y(pt2))
+ eq++;
+ }
+ CU_ASSERT_EQUAL(eq, 1000);
+ lwmpoint_free(mpt2);
+ pt = NULL;
+ pt2 = NULL;
+
+
+ /* 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);
+ // 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);
lwmpoint_free(mpt);
+ pt = NULL;
- mpt = lwgeom_to_points(geom, 0);
+ mpt = lwgeom_to_points(geom, 0, 0);
CU_ASSERT_EQUAL(mpt, NULL);
lwmpoint_free(mpt);
/* MULTIPOLYGON */
geom = lwgeom_from_wkt("MULTIPOLYGON(((10 0,0 10,10 20,20 10,10 0)),((0 0,5 0,5 5,0 5,0 0)))", LW_PARSER_CHECK_NONE);
- mpt = lwgeom_to_points(geom, 1000);
+ mpt = lwgeom_to_points(geom, 1000, 0);
CU_ASSERT_EQUAL(mpt->ngeoms,1000);
lwmpoint_free(mpt);
- mpt = lwgeom_to_points(geom, 1);
+ mpt = lwgeom_to_points(geom, 1, 0);
CU_ASSERT_EQUAL(mpt->ngeoms,1);
lwmpoint_free(mpt);
/*
* Point density functions
*/
-extern LWMPOINT *lwpoly_to_points(const LWPOLY *poly, uint32_t npoints);
-extern LWMPOINT *lwmpoly_to_points(const LWMPOLY *mpoly, uint32_t npoints);
-extern LWMPOINT *lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints);
+extern LWMPOINT *lwpoly_to_points(const LWPOLY *poly, uint32_t npoints, int32_t seed);
+extern LWMPOINT *lwmpoly_to_points(const LWMPOLY *mpoly, uint32_t npoints, int32_t seed);
+extern LWMPOINT *lwgeom_to_points(const LWGEOM *lwgeom, uint32_t npoints, int32_t seed);
/*
* Geometric median
#include "liblwgeom.h"
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
+#include "lwrandom.h"
#include <stdarg.h>
#include <stdlib.h>
-#include <time.h>
LWTIN* lwtin_from_geos(const GEOSGeometry* geom, uint8_t want3d);
}
LWMPOINT*
-lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints)
+lwpoly_to_points(const LWPOLY* lwpoly, uint32_t npoints, int32_t seed)
{
double area, bbox_area, bbox_width, bbox_height;
GBOX bbox;
* 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;
+ if (sample_sqrt == 0)
+ sample_sqrt = 1;
/* Calculate the grids we're going to randomize within */
if (bbox_width > bbox_height)
/* Get an empty multi-point ready to return */
mpt = lwmpoint_construct_empty(srid, 0, 0);
- /* Init random number generator */
- srand(time(NULL));
+ /* Initiate random number generator.
+ * Repeatable numbers are generated with seed values >= 1.
+ * When seed is zero and has not previously been set, it is based on
+ * 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 */
}
}
- /* shuffle */
+ /* Fisher-Yates shuffle */
n = sample_height * sample_width;
if (n > 1)
{
- for (i = 0; i < n - 1; ++i)
+ for (i = n - 1; i > 0; i--)
{
- size_t rnd = (size_t)rand();
- size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
+ 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);
int contains = 0;
double y = bbox.ymin + cells[2 * i] * sample_cell_size;
double x = bbox.xmin + cells[2 * i + 1] * sample_cell_size;
- x += rand() * sample_cell_size / RAND_MAX;
- y += rand() * sample_cell_size / RAND_MAX;
+ 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);
/* Allocate points to sub-geometries by area, then call lwgeom_poly_to_points and bundle up final result in a single
* multipoint. */
LWMPOINT*
-lwmpoly_to_points(const LWMPOLY* lwmpoly, uint32_t npoints)
+lwmpoly_to_points(const LWMPOLY* lwmpoly, uint32_t npoints, int32_t seed)
{
const LWGEOM* lwgeom = (LWGEOM*)lwmpoly;
double area;
int sub_npoints = lround(npoints * sub_area / area);
if (sub_npoints > 0)
{
- LWMPOINT* sub_mpt = lwpoly_to_points(lwmpoly->geoms[i], sub_npoints);
+ LWMPOINT* sub_mpt = lwpoly_to_points(lwmpoly->geoms[i], sub_npoints, seed);
if (!mpt)
mpt = sub_mpt;
else
}
LWMPOINT*
-lwgeom_to_points(const LWGEOM* lwgeom, uint32_t npoints)
+lwgeom_to_points(const LWGEOM* lwgeom, uint32_t npoints, int32_t seed)
{
switch (lwgeom_get_type(lwgeom))
{
case MULTIPOLYGONTYPE:
- return lwmpoly_to_points((LWMPOLY*)lwgeom, npoints);
+ return lwmpoly_to_points((LWMPOLY*)lwgeom, npoints, seed);
case POLYGONTYPE:
- return lwpoly_to_points((LWPOLY*)lwgeom, npoints);
+ return lwpoly_to_points((LWPOLY*)lwgeom, npoints, seed);
default:
lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(lwgeom_get_type(lwgeom)));
return NULL;
--- /dev/null
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * Copyright 2019 Mike Taves
+ *
+ **********************************************************************/
+
+#include "lwrandom.h"
+
+#include <stdlib.h>
+#include <stdint.h>
+#include <stddef.h>
+#include <time.h>
+
+#ifdef _WIN32
+#include <process.h>
+#define getpid _getpid
+#else
+#include <unistd.h>
+#endif
+
+static unsigned char _lwrandom_seed_set = 0;
+static int32_t _lwrandom_seed[3] = {0x330e, 0xabcd, 0x1234};
+
+/*
+ * Set seed for a random number generator.
+ * Repeatable numbers are generated with seed values >= 1.
+ * When seed is zero and has not previously been set, it is based on
+ * Unix time (seconds) and process ID. */
+void
+lwrandom_set_seed(int32_t seed)
+{
+ if (seed == 0)
+ {
+ if (_lwrandom_seed_set == 0)
+ seed = (unsigned int)time(NULL) + (unsigned int)getpid() - 0xbadd;
+ else
+ return;
+ }
+ /* s1 value between 1 and 2147483562 */
+ _lwrandom_seed[1] = (((int64_t)seed + 0xfeed) % 2147483562) + 1;
+ /* s2 value between 1 and 2147483398 */
+ _lwrandom_seed[2] = ((((int64_t)seed + 0xdefeb) << 5) % 2147483398) + 1;
+ _lwrandom_seed_set = 1;
+}
+
+/* for low-level external debugging */
+void
+_lwrandom_set_seeds(int32_t s1, int32_t s2)
+{
+ /* _lwrandom_seed[0] not used */
+ _lwrandom_seed[1] = s1;
+ _lwrandom_seed[2] = s2;
+ _lwrandom_seed_set = 1;
+}
+int32_t
+_lwrandom_get_seed(size_t idx)
+{
+ return _lwrandom_seed[idx];
+}
+
+/*
+ * Generate a random floating-point value.
+ * Values are uniformly distributed between 0 and 1.
+ *
+ * Authors:
+ * Pierre L'Ecuyer (1988), see source code in Figure 3.
+ * C version by John Burkardt, modified by Mike Taves.
+ *
+ * Reference:
+ * Pierre L'Ecuyer,
+ * Efficient and Portable Combined Random Number Generators,
+ * Communications of the ACM, Volume 31, Number 6, June 1988,
+ * pages 742-751. doi:10.1145/62959.62969
+ */
+double
+lwrandom_uniform(void)
+{
+ double value;
+ int32_t k;
+ int32_t z;
+ int32_t *s1 = &_lwrandom_seed[1];
+ int32_t *s2 = &_lwrandom_seed[2];
+
+ k = *s1 / 53668;
+ *s1 = 40014 * (*s1 - k * 53668) - k * 12211;
+ if (*s1 < 0)
+ *s1 += 2147483563;
+
+ k = *s2 / 52774;
+ *s2 = 40692 * (*s2 - k * 52774) - k * 3791;
+ if (*s2 < 0)
+ *s2 += 2147483399;
+
+ z = *s1 - *s2;
+ if (z < 1)
+ z += 2147483562;
+
+ value = (double)(z) / 2147483563.0;
+
+ return value;
+}
--- /dev/null
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * PostGIS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * PostGIS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ **********************************************************************
+ *
+ * Copyright 2019 Mike Taves
+ *
+ **********************************************************************/
+
+#include <stddef.h>
+#include <stdint.h>
+
+void lwrandom_set_seed(int32_t seed);
+double lwrandom_uniform(void);
+
+/* for low-level external debugging */
+void _lwrandom_set_seeds(int32_t s1, int32_t s2);
+int32_t _lwrandom_get_seed(size_t idx);
LWGEOM *lwgeom_input;
LWGEOM *lwgeom_result;
int32 npoints;
+ int32 seed = 0;
gser_input = PG_GETARG_GSERIALIZED_P(0);
- npoints = DatumGetInt32(DirectFunctionCall1(numeric_int4, PG_GETARG_DATUM(1)));
+ npoints = PG_GETARG_INT32(1);
if (npoints < 0)
PG_RETURN_NULL();
+ if (PG_NARGS() > 2 && ! PG_ARGISNULL(2))
+ {
+ seed = PG_GETARG_INT32(2);
+ if (seed < 1)
+ {
+ lwpgerror("ST_GeneratePoints: seed must be greater than zero");
+ PG_RETURN_NULL();
+ }
+ }
+
/* Types get checked in the code, we'll keep things small out there */
lwgeom_input = lwgeom_from_gserialized(gser_input);
- lwgeom_result = (LWGEOM*)lwgeom_to_points(lwgeom_input, npoints);
+ lwgeom_result = (LWGEOM*)lwgeom_to_points(lwgeom_input, npoints, seed);
lwgeom_free(lwgeom_input);
PG_FREE_IF_COPY(gser_input, 0);
COST 1; -- reset cost, see #3675
-- Availability: 2.3.0
-CREATE OR REPLACE FUNCTION ST_GeneratePoints(area geometry, npoints numeric)
+CREATE OR REPLACE FUNCTION ST_GeneratePoints(area geometry, npoints integer)
RETURNS geometry
AS 'MODULE_PATHNAME','ST_GeneratePoints'
- LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL
+ LANGUAGE 'c' VOLATILE STRICT _PARALLEL
COST 1; -- reset cost, see #3675
+-- Availability: 3.0.0
+CREATE OR REPLACE FUNCTION ST_GeneratePoints(area geometry, npoints integer, seed integer)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME','ST_GeneratePoints'
+ LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL
+ COST 1;
+
-- PostGIS equivalent function: convexhull(geometry)
CREATE OR REPLACE FUNCTION ST_ConvexHull(geometry)
RETURNS geometry
#if POSTGIS_PGSQL_VERSION >= 96
#define _PARALLEL PARALLEL SAFE
+#define _PARALLEL_RESTRICTED PARALLEL RESTRICTED
#else
#define _PARALLEL
+#define _PARALLEL_RESTRICTED
#endif
/*
ANALYZE bug_4144_table;
DROP TABLE IF EXISTS bug_4144_table;
+-- #4299
+SELECT '#4299', ST_Disjoint(ST_GeneratePoints(g, 1000), ST_GeneratePoints(g, 1000))
+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
+FROM (SELECT 'POLYGON((0 0,1 0,1 1,0 1,0 0))'::geometry AS g) AS f;
+
+
-- Clean up
DELETE FROM spatial_ref_sys;
ERROR: lwgeom_pointonsurface: GEOS Error: TopologyException: Input geom 1 is invalid: Self-intersection
#4081|f|t
NOTICE: table "bug_4144_table" does not exist, skipping
+#4299|t
+#4304|t|t|t|t
#4176|t
projects / postgis / commitdiff