[postgis-tickets] r14492 - #3399, ST_GeneratePoints
Paul Ramsey
pramsey at cleverelephant.ca
Thu Dec 17 10:12:24 PST 2015
Author: pramsey
Date: 2015-12-17 10:12:23 -0800 (Thu, 17 Dec 2015)
New Revision: 14492
Modified:
trunk/NEWS
trunk/doc/reference_processing.xml
trunk/liblwgeom/cunit/cu_algorithm.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/lwgeom_geos.c
trunk/postgis/lwgeom_geos.c
trunk/postgis/postgis.sql.in
Log:
#3399, ST_GeneratePoints
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/NEWS 2015-12-17 18:12:23 UTC (rev 14492)
@@ -10,6 +10,7 @@
- TopoGeom_addElement, TopoGeom_remElement (Sandro Santilli)
- populate_topology_layer (Sandro Santilli)
- #2259 ST_Voronoi (Dan Baston)
+ - #3339 ST_GeneratePoints (Paul Ramsey)
PostGIS 2.2.0
2015/10/07
Modified: trunk/doc/reference_processing.xml
===================================================================
--- trunk/doc/reference_processing.xml 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/doc/reference_processing.xml 2015-12-17 18:12:23 UTC (rev 14492)
@@ -3653,4 +3653,43 @@
</para>
</refsection>
</refentry>
+
+ <refentry id="ST_GeneratePoints">
+ <refnamediv>
+ <refname>ST_GeneratePoints</refname>
+
+ <refpurpose>Converts a polygon or multi-polygon into a multi-point composed of randomly location points within the original areas.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_GeneratePoints</function></funcdef>
+ <paramdef>
+ <parameter>g</parameter>
+ <type>geometry</type>
+ </paramdef>
+ <paramdef>
+ <parameter>npoints</parameter>
+ <type>numeric</type>
+ </paramdef>
+ </funcprototype>
+
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>
+ ST_GeneratePoints generates pseudo-random points until the requested number are
+ found within the input area.
+ </para>
+
+ <para>Availability: 2.3.0</para>
+ </refsection>
+
+ </refentry>
+
+
</sect1>
Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/liblwgeom/cunit/cu_algorithm.c 2015-12-17 18:12:23 UTC (rev 14492)
@@ -982,6 +982,47 @@
lwfree(ewkt);
}
+
+
+static void test_point_density(void)
+{
+ LWGEOM *geom;
+ LWMPOINT *mpt;
+ // 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);
+ CU_ASSERT_EQUAL(mpt->ngeoms,100);
+ // ewkt = lwgeom_to_ewkt((LWGEOM*)mpt);
+ // printf("%s\n", ewkt);
+ // lwfree(ewkt);
+ lwmpoint_free(mpt);
+
+ mpt = lwgeom_to_points(geom, 1);
+ CU_ASSERT_EQUAL(mpt->ngeoms,1);
+ lwmpoint_free(mpt);
+
+ mpt = lwgeom_to_points(geom, 0);
+ CU_ASSERT_EQUAL(mpt, NULL);
+
+ lwgeom_free(geom);
+
+ /* 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);
+ CU_ASSERT_EQUAL(mpt->ngeoms,1000);
+ lwmpoint_free(mpt);
+
+ mpt = lwgeom_to_points(geom, 1);
+ CU_ASSERT_EQUAL(mpt->ngeoms,1);
+ lwmpoint_free(mpt);
+
+ lwgeom_free(geom);
+}
+
+
/*
** Used by test harness to register the tests in this file.
*/
@@ -1007,4 +1048,5 @@
PG_ADD_TEST(suite,test_isclosed);
PG_ADD_TEST(suite,test_lwgeom_simplify);
PG_ADD_TEST(suite,test_lw_arc_center);
+ PG_ADD_TEST(suite,test_point_density);
}
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/liblwgeom/liblwgeom.h.in 2015-12-17 18:12:23 UTC (rev 14492)
@@ -1440,6 +1440,13 @@
extern LWPOLY *lwpoly_segmentize2d(LWPOLY *line, double dist);
extern LWCOLLECTION *lwcollection_segmentize2d(LWCOLLECTION *coll, double dist);
+/*
+ * Point density functions
+ */
+extern LWMPOINT *lwpoly_to_points(const LWPOLY *poly, int npoints);
+extern LWMPOINT *lwmpoly_to_points(const LWMPOLY *mpoly, int npoints);
+extern LWMPOINT *lwgeom_to_points(const LWGEOM *lwgeom, int npoints);
+
/**
* Calculate the GeoHash (http://geohash.org) string for a geometry. Caller must free.
*/
Modified: trunk/liblwgeom/lwgeom_geos.c
===================================================================
--- trunk/liblwgeom/lwgeom_geos.c 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/liblwgeom/lwgeom_geos.c 2015-12-17 18:12:23 UTC (rev 14492)
@@ -29,6 +29,7 @@
#include "lwgeom_log.h"
#include <stdlib.h>
+#include <time.h>
LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d);
@@ -1546,6 +1547,265 @@
return lwgeom_result;
}
+
+static void shuffle(void *array, size_t n, size_t size) {
+ char tmp[size];
+ char *arr = array;
+ size_t stride = size;
+
+ if (n > 1) {
+ size_t i;
+ for (i = 0; i < n - 1; ++i) {
+ size_t rnd = (size_t) rand();
+ size_t j = i + rnd / (RAND_MAX / (n - i) + 1);
+
+ memcpy(tmp, arr + j * stride, size);
+ memcpy(arr + j * stride, arr + i * stride, size);
+ memcpy(arr + i * stride, tmp, size);
+ }
+ }
+}
+
+LWMPOINT*
+lwpoly_to_points(const LWPOLY *lwpoly, int npoints)
+{
+ double area, bbox_area, bbox_width, bbox_height;
+ GBOX bbox;
+ const LWGEOM *lwgeom = (LWGEOM*)lwpoly;
+ int sample_npoints, sample_sqrt, sample_width, sample_height;
+ double sample_cell_size;
+ int i, j;
+ int iterations = 0;
+ int npoints_generated = 0;
+ int npoints_tested = 0;
+ GEOSGeometry *g;
+ const GEOSPreparedGeometry *gprep;
+ GEOSGeometry *gpt;
+ GEOSCoordSequence *gseq;
+ LWMPOINT *mpt;
+ int srid = lwgeom_get_srid(lwgeom);
+ int done = 0;
+ int *cells;
+
+ if (lwgeom_get_type(lwgeom) != POLYGONTYPE)
+ {
+ lwerror("%s: only polygons supported", __func__);
+ return NULL;
+ }
+
+ if (npoints == 0 || lwgeom_is_empty(lwgeom))
+ {
+ return NULL;
+ // return lwmpoint_construct_empty(lwgeom_get_srid(poly), lwgeom_has_z(poly), lwgeom_has_m(poly));
+ }
+
+ if (!lwpoly->bbox)
+ {
+ lwgeom_calculate_gbox(lwgeom, &bbox);
+ }
+ else
+ {
+ bbox = *(lwpoly->bbox);
+ }
+ area = lwpoly_area(lwpoly);
+ bbox_width = bbox.xmax - bbox.xmin;
+ bbox_height = bbox.ymax - bbox.ymin;
+ bbox_area = bbox_width * bbox_height;
+
+ if (area == 0.0 || bbox_area == 0.0)
+ {
+ lwerror("%s: zero area input polygon, TBD", __func__);
+ return NULL;
+ }
+
+ /* Gross up our test set a bit to increase odds of getting */
+ /* coverage in one pass */
+ sample_npoints = npoints * bbox_area / area;
+
+ /* 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;
+
+ /* Calculate the grids we're going to randomize within */
+ if (bbox_width > bbox_height)
+ {
+ sample_width = sample_sqrt;
+ sample_height = ceil((double)sample_npoints / (double)sample_width);
+ sample_cell_size = bbox_width / sample_width;
+ }
+ else
+ {
+ sample_height = sample_sqrt;
+ sample_width = ceil((double)sample_npoints / (double)sample_height);
+ sample_cell_size = bbox_height / sample_height;
+ }
+
+ /* Prepare the polygon for fast true/false testing */
+ initGEOS(lwnotice, lwgeom_geos_error);
+ g = (GEOSGeometry *)LWGEOM2GEOS(lwgeom, 0);
+ if (!g)
+ {
+ lwerror("%s: Geometry could not be converted to GEOS: %s", __func__, lwgeom_geos_errmsg);
+ return NULL;
+ }
+ gprep = GEOSPrepare(g);
+
+ /* Get an empty multi-point ready to return */
+ mpt = lwmpoint_construct_empty(srid, 0, 0);
+
+ /* Init random number generator */
+ srand(time(NULL));
+
+ /* 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;
+ }
+ }
+ shuffle(cells, sample_height*sample_width, 2*sizeof(int));
+
+ /* Start testing points */
+ while (npoints_generated < npoints)
+ {
+ iterations++;
+ for (i = 0; i < sample_width*sample_height; 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;
+ x += rand() * sample_cell_size / RAND_MAX;
+ y += rand() * sample_cell_size / RAND_MAX;
+ if (x >= bbox.xmax || y >= bbox.ymax)
+ continue;
+
+ gseq = GEOSCoordSeq_create(1, 2);
+ GEOSCoordSeq_setX(gseq, 0, x);
+ GEOSCoordSeq_setY(gseq, 0, y);
+ gpt = GEOSGeom_createPoint(gseq);
+
+ contains = GEOSPreparedIntersects(gprep, gpt);
+
+ GEOSGeom_destroy(gpt);
+
+ if (contains == 2)
+ {
+ GEOSPreparedGeom_destroy(gprep);
+ GEOSGeom_destroy(g);
+ lwerror("%s: GEOS exception on PreparedContains: %s", __func__, lwgeom_geos_errmsg);
+ return NULL;
+ }
+ if (contains)
+ {
+ npoints_generated++;
+ mpt = lwmpoint_add_lwpoint(mpt, lwpoint_make2d(srid, x, y));
+ if (npoints_generated == npoints)
+ {
+ done = 1;
+ break;
+ }
+ }
+
+ /* Short-circuit check for ctrl-c occasionally */
+ npoints_tested++;
+ if (npoints_tested % 10000 == 0)
+ {
+ LW_ON_INTERRUPT(GEOSPreparedGeom_destroy(gprep); GEOSGeom_destroy(g); return NULL);
+ }
+
+ if (done) break;
+ }
+ if (done || iterations > 100) break;
+ }
+
+ GEOSPreparedGeom_destroy(gprep);
+ GEOSGeom_destroy(g);
+ lwfree(cells);
+
+ return mpt;
+}
+
+
+/*
+* 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, int npoints)
+{
+ const LWGEOM *lwgeom = (LWGEOM*)lwmpoly;
+ double area;
+ int i;
+ LWMPOINT *mpt = NULL;
+
+ if (lwgeom_get_type(lwgeom) != MULTIPOLYGONTYPE)
+ {
+ lwerror("%s: only multipolygons supported", __func__);
+ return NULL;
+ }
+ if (npoints == 0 || lwgeom_is_empty(lwgeom))
+ {
+ return NULL;
+ }
+
+ area = lwgeom_area(lwgeom);
+
+ for (i = 0; i < lwmpoly->ngeoms; i++)
+ {
+ double sub_area = lwpoly_area(lwmpoly->geoms[i]);
+ int sub_npoints = lround(npoints * sub_area / area);
+ if(sub_npoints > 0)
+ {
+ LWMPOINT *sub_mpt = lwpoly_to_points(lwmpoly->geoms[i], sub_npoints);
+ if (!mpt)
+ {
+ mpt = sub_mpt;
+ }
+ else
+ {
+ int j;
+ for (j = 0; j < sub_mpt->ngeoms; j++)
+ {
+ mpt = lwmpoint_add_lwpoint(mpt, sub_mpt->geoms[j]);
+ }
+ /* Just free the shell, leave the underlying lwpoints alone, as they
+ are now owed by the returning multipoint */
+ lwgeom_release((LWGEOM*)sub_mpt);
+ }
+ }
+ }
+
+ return mpt;
+}
+
+
+LWMPOINT*
+lwgeom_to_points(const LWGEOM *lwgeom, int npoints)
+{
+ switch(lwgeom_get_type(lwgeom))
+ {
+ case MULTIPOLYGONTYPE:
+ return lwmpoly_to_points((LWMPOLY*)lwgeom, npoints);
+ case POLYGONTYPE:
+ return lwpoly_to_points((LWPOLY*)lwgeom, npoints);
+ default:
+ lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(lwgeom_get_type(lwgeom)));
+ return NULL;
+ }
+}
+
+
+
+
LWTIN *lwtin_from_geos(const GEOSGeometry *geom, int want3d) {
int type = GEOSGeomTypeId(geom);
int hasZ;
Modified: trunk/postgis/lwgeom_geos.c
===================================================================
--- trunk/postgis/lwgeom_geos.c 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/postgis/lwgeom_geos.c 2015-12-17 18:12:23 UTC (rev 14492)
@@ -33,6 +33,7 @@
#include "utils/array.h"
#include "utils/builtins.h"
#include "utils/lsyscache.h"
+#include "utils/numeric.h"
#if POSTGIS_PGSQL_VERSION >= 93
#include "access/htup_details.h"
@@ -917,7 +918,48 @@
PG_RETURN_POINTER(result);
}
+
+
/*
+* Generate a field of random points within the area of a
+* polygon or multipolygon. Throws an error for other geometry
+* types.
+*/
+Datum ST_GeneratePoints(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_GeneratePoints);
+Datum ST_GeneratePoints(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gser_input;
+ GSERIALIZED *gser_result;
+ LWGEOM *lwgeom_input;
+ LWGEOM *lwgeom_result;
+ int32 npoints;
+
+ gser_input = PG_GETARG_GSERIALIZED_P(0);
+ npoints = DatumGetInt32(DirectFunctionCall1(numeric_int4, PG_GETARG_DATUM(1)));
+
+ /* Smartasses get nothing back */
+ if (npoints < 0)
+ 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_free(lwgeom_input);
+ PG_FREE_IF_COPY(gser_input, 0);
+
+ /* Return null as null */
+ if (!lwgeom_result)
+ PG_RETURN_NULL();
+
+ /* Serialize and return */
+ gser_result = gserialized_from_lwgeom(lwgeom_result, 0);
+ lwgeom_free(lwgeom_result);
+ PG_RETURN_POINTER(gser_result);
+}
+
+
+/*
* Compute at offset curve to a line
*/
Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2015-12-17 07:30:43 UTC (rev 14491)
+++ trunk/postgis/postgis.sql.in 2015-12-17 18:12:23 UTC (rev 14492)
@@ -3197,6 +3197,13 @@
LANGUAGE 'c' IMMUTABLE STRICT
COST 100;
+-- Availability: 2.3.0
+CREATE OR REPLACE FUNCTION ST_GeneratePoints(area geometry, npoints numeric)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME','ST_GeneratePoints'
+ LANGUAGE 'c' IMMUTABLE STRICT
+ COST 400;
+
-- PostGIS equivalent function: convexhull(geometry)
CREATE OR REPLACE FUNCTION ST_ConvexHull(geometry)
RETURNS geometry
More information about the postgis-tickets
mailing list