[postgis-tickets] r14467 - #2259, ST_Voronoi
Daniel Baston
dbaston at gmail.com
Tue Dec 1 17:10:34 PST 2015
Author: dbaston
Date: 2015-12-01 17:10:33 -0800 (Tue, 01 Dec 2015)
New Revision: 14467
Added:
trunk/regress/voronoi.sql
trunk/regress/voronoi_expected
Modified:
trunk/doc/reference_processing.xml
trunk/liblwgeom/cunit/cu_triangulate.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/lwgeom_geos.c
trunk/postgis/lwgeom_geos.c
trunk/postgis/postgis.sql.in
trunk/regress/Makefile.in
Log:
#2259, ST_Voronoi
Modified: trunk/doc/reference_processing.xml
===================================================================
--- trunk/doc/reference_processing.xml 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/doc/reference_processing.xml 2015-12-02 01:10:33 UTC (rev 14467)
@@ -3577,4 +3577,76 @@
</para>
</refsection>
</refentry>
+
+ <refentry id="ST_Voronoi">
+ <refnamediv>
+ <refname>ST_Voronoi</refname>
+
+ <refpurpose>Computes a Voronoi diagram from the vertices of a geometry.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_Voronoi</function></funcdef>
+ <paramdef>
+ <parameter>g1</parameter>
+ <type>geometry</type>
+ </paramdef>
+ <paramdef choice="opt">
+ <parameter>clip</parameter>
+ <type>geometry</type>
+ </paramdef>
+ <paramdef choice="opt">
+ <parameter>tolerance</parameter>
+ <type>float8</type>
+ </paramdef>
+ <paramdef choice="opt">
+ <parameter>return_polygons</parameter>
+ <type>boolean</type>
+ </paramdef>
+ </funcprototype>
+
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>
+ ST_Voronoi computes a two-dimensional Voronoi diagram from the vertices of
+ the supplied geometry. By default, the result will be a GeometryCollection of
+ Polygons that covers an envelope larger than the extent of the input vertices.
+ </para>
+
+ <para>
+ Optional parameters:
+ <itemizedlist>
+ <listitem>
+ <para>'clip' : If a geometry is supplied as the "clip" parameter, the diagram will be extended to cover the
+ envelope of the "clip" geometry, unless that envelope is smaller than the default envelope.
+ (default = NULL)</para>
+ </listitem>
+ <listitem>
+ <para>'tolerance' : The distance within which vertices will be considered equivalent. Robustness of the algorithm can be improved by supplying a nonzero tolerance distance. (default = 0.0)</para>
+ </listitem>
+ <listitem>
+ <para>'return_polygons' : if true, the result of ST_Voronoi will be a GeometryCollection of Polygons. If false, the result will be a MultiLineString. (default = true)</para>
+ </listitem>
+ </itemizedlist>
+ </para>
+
+ <para>Availability: 2.3.0 - requires GEOS >= 3.5.0.</para>
+ </refsection>
+
+ <!-- Optionally add a "See Also" section -->
+ <refsection>
+ <title>See Also</title>
+
+ <para>
+ <xref linkend="ST_DelaunayTriangles" />,
+ <xref linkend="ST_Collect" />
+ </para>
+ </refsection>
+ </refentry>
</sect1>
Modified: trunk/liblwgeom/cunit/cu_triangulate.c
===================================================================
--- trunk/liblwgeom/cunit/cu_triangulate.c 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/liblwgeom/cunit/cu_triangulate.c 2015-12-02 01:10:33 UTC (rev 14467)
@@ -3,6 +3,7 @@
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
*
+ * Copyright (C) 2015 Daniel Baston <dbaston at gmail.com>
* Copyright (C) 2012 Sandro Santilli <strk at keybit.net>
*
* This is free software; you can redistribute and/or modify it under
@@ -13,7 +14,6 @@
#include "CUnit/Basic.h"
#include "cu_tester.h"
-#include "liblwgeom.h"
#include "liblwgeom_internal.h"
static void test_lwgeom_delaunay_triangulation(void)
@@ -43,10 +43,70 @@
CU_ASSERT_STRING_EQUAL(wkt, exp_wkt);
lwfree(wkt);
-#endif /* POSTGIS_GEOS_VERSION >= 33 */
+#endif /* POSTGIS_GEOS_VERSION >= 34 */
}
+static void test_lwgeom_voronoi_diagram(void)
+{
+#if POSTGIS_GEOS_VERSION >= 35
+ LWGEOM* in = lwgeom_from_wkt("MULTIPOINT(4 4, 5 5, 6 6)", LW_PARSER_CHECK_NONE);
+ LWGEOM* out_boundaries = lwgeom_voronoi_diagram(in, NULL, 0, 0);
+ LWGEOM* out_lines = lwgeom_voronoi_diagram(in, NULL, 0, 1);
+
+ /* For boundaries we get a generic LWCOLLECTION */
+ CU_ASSERT_EQUAL(COLLECTIONTYPE, lwgeom_get_type(out_boundaries));
+ /* For lines we get a MULTILINETYPE */
+ CU_ASSERT_EQUAL(MULTILINETYPE, lwgeom_get_type(out_lines));
+
+ lwgeom_free(in);
+ lwgeom_free(out_boundaries);
+ lwgeom_free(out_lines);
+#endif /* POSTGIS_GEOS_VERSION >= 35 */
+}
+
+static void test_lwgeom_voronoi_diagram_custom_envelope(void)
+{
+#if POSTGIS_GEOS_VERSION >= 35
+ LWGEOM* in = lwgeom_from_wkt("MULTIPOINT(4 4, 5 5, 6 6)", LW_PARSER_CHECK_NONE);
+ LWGEOM* for_extent = lwgeom_from_wkt("LINESTRING (-10 -10, 10 10)", LW_PARSER_CHECK_NONE);
+ const GBOX* clipping_extent = lwgeom_get_bbox(for_extent);
+
+ LWGEOM* out = lwgeom_voronoi_diagram(in, clipping_extent, 0, 0);
+ const GBOX* output_extent = lwgeom_get_bbox(out);
+
+ CU_ASSERT_TRUE(gbox_same(clipping_extent, output_extent));
+
+ lwgeom_free(in);
+ lwgeom_free(for_extent);
+ lwgeom_free(out);
+#endif /* POSTGIS_GEOS_VERSION >= 35 */
+}
+
+static void assert_empty_diagram(char* wkt, double tolerance)
+{
+ LWGEOM* in = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_NONE);
+ LWGEOM* out = lwgeom_voronoi_diagram(in, NULL, 0, 0);
+
+ CU_ASSERT_TRUE(lwgeom_is_collection(out));
+ CU_ASSERT_EQUAL(COLLECTIONTYPE, lwgeom_get_type(out));
+
+ lwgeom_free(in);
+ lwgeom_free(out);
+}
+
+static void test_lwgeom_voronoi_diagram_expected_empty(void)
+{
+#if POSTGIS_GEOS_VERSION >= 35
+ assert_empty_diagram("POLYGON EMPTY", 0);
+ assert_empty_diagram("POINT (1 2)", 0);
+
+ /* This one produces an empty diagram because our two unqiue points
+ * collapse onto one after considering the tolerance. */
+ assert_empty_diagram("MULTIPOINT (0 0, 0 0.00001)", 0.001);
+#endif /* POSTGIS_GEOS_VERSION >= 35 */
+}
+
/*
** Used by test harness to register the tests in this file.
*/
@@ -55,4 +115,7 @@
{
CU_pSuite suite = CU_add_suite("triangulate", NULL, NULL);
PG_ADD_TEST(suite, test_lwgeom_delaunay_triangulation);
+ PG_ADD_TEST(suite, test_lwgeom_voronoi_diagram);
+ PG_ADD_TEST(suite, test_lwgeom_voronoi_diagram_expected_empty);
+ PG_ADD_TEST(suite, test_lwgeom_voronoi_diagram_custom_envelope);
}
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/liblwgeom/liblwgeom.h.in 2015-12-02 01:10:33 UTC (rev 14467)
@@ -2238,11 +2238,24 @@
* triangulation on them.
*
* @param geom the input geometry
- * @param tolerance an optional snapping tolerance for improved tolerance
+ * @param tolerance an optional snapping tolerance for improved robustness
* @param edgeOnly if non-zero the result will be a MULTILINESTRING,
* otherwise it'll be a COLLECTION of polygons.
*/
LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *geom, double tolerance, int edgeOnly);
+/**
+ * Take vertices of a geometry and build the Voronoi diagram
+ *
+ * @param g the input geometry
+ * @param env an optional envelope for clipping the results
+ * @param tolerance an optional snapping tolerance for improved robustness
+ * @param output_edges if non-zero the result will be a MULTILINESTRING,
+ * otherwise it'll be a COLLECTION of polygons.
+ *
+ * Requires GEOS-3.5.0 or higher
+ */
+LWGEOM* lwgeom_voronoi_diagram(const LWGEOM* g, const GBOX* env, double tolerance, int output_edges);
+
#endif /* !defined _LIBLWGEOM_H */
Modified: trunk/liblwgeom/lwgeom_geos.c
===================================================================
--- trunk/liblwgeom/lwgeom_geos.c 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/liblwgeom/lwgeom_geos.c 2015-12-02 01:10:33 UTC (rev 14467)
@@ -1657,3 +1657,101 @@
#endif /* POSTGIS_GEOS_VERSION < 34 */
}
+
+static
+GEOSCoordSequence* lwgeom_get_geos_coordseq_2d(const LWGEOM* g, uint32_t num_points)
+{
+ uint32_t i = 0;
+ uint8_t num_dims = 2;
+ LWPOINTITERATOR* it;
+ GEOSCoordSequence* coords;
+ POINT4D tmp;
+
+ coords = GEOSCoordSeq_create(num_points, num_dims);
+ if (!coords)
+ return NULL;
+
+ it = lwpointiterator_create(g);
+ while(lwpointiterator_next(it, &tmp))
+ {
+ if(i >= num_points)
+ {
+ lwerror("Incorrect num_points provided to lwgeom_get_geos_coordseq_2d");
+ GEOSCoordSeq_destroy(coords);
+ lwpointiterator_destroy(it);
+ return NULL;
+ }
+
+ if(!GEOSCoordSeq_setX(coords, i, tmp.x) || !GEOSCoordSeq_setY(coords, i, tmp.y))
+ {
+ GEOSCoordSeq_destroy(coords);
+ lwpointiterator_destroy(it);
+ return NULL;
+ }
+ i++;
+ }
+ lwpointiterator_destroy(it);
+
+ return coords;
+}
+
+LWGEOM* lwgeom_voronoi_diagram(const LWGEOM* g, const GBOX* env, double tolerance, int output_edges) {
+#if POSTGIS_GEOS_VERSION < 35
+ lwerror("lwgeom_voronoi_diagram: GEOS 3.5 or higher required");
+ return NULL;
+#else
+ uint32_t num_points = lwgeom_count_vertices(g);
+ LWGEOM *lwgeom_result;
+ char is_3d = LW_FALSE;
+ int srid = lwgeom_get_srid(g);
+ GEOSCoordSequence* coords;
+ GEOSGeometry* geos_geom;
+ GEOSGeometry* geos_env = NULL;
+ GEOSGeometry* geos_result;
+
+ if (num_points < 2)
+ {
+ LWCOLLECTION* empty = lwcollection_construct_empty(COLLECTIONTYPE, lwgeom_get_srid(g), 0, 0);
+ return lwcollection_as_lwgeom(empty);
+ }
+
+ initGEOS(lwnotice, lwgeom_geos_error);
+
+ /* Instead of using the standard LWGEOM2GEOS transformer, we read the vertices of the
+ * LWGEOM directly and put them into a single GEOS CoordinateSeq that can be used to
+ * define a LineString. This allows us to process geometry types that may not be
+ * supported by GEOS, and reduces the memory requirements in cases of many geometries
+ * with few points (such as LWMPOINT).
+ */
+ coords = lwgeom_get_geos_coordseq_2d(g, num_points);
+ if (!coords)
+ return NULL;
+
+ geos_geom = GEOSGeom_createLineString(coords);
+ if (!geos_geom)
+ {
+ GEOSCoordSeq_destroy(coords);
+ return NULL;
+ }
+
+ if (env)
+ geos_env = GBOX2GEOS(env);
+
+ geos_result = GEOSVoronoiDiagram(geos_geom, geos_env, tolerance, output_edges);
+
+ GEOSGeom_destroy(geos_geom);
+ if (env)
+ GEOSGeom_destroy(geos_env);
+
+ if (!geos_result)
+ return NULL;
+
+ lwgeom_result = GEOS2LWGEOM(geos_result, is_3d);
+ GEOSGeom_destroy(geos_result);
+
+ lwgeom_set_srid(lwgeom_result, srid);
+
+ return lwgeom_result;
+#endif /* POSTGIS_GEOS_VERSION < 35 */
+}
+
Modified: trunk/postgis/lwgeom_geos.c
===================================================================
--- trunk/postgis/lwgeom_geos.c 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/postgis/lwgeom_geos.c 2015-12-02 01:10:33 UTC (rev 14467)
@@ -3441,3 +3441,102 @@
PG_RETURN_POINTER(out);
}
+/******************************************
+ *
+ * ST_Voronoi
+ *
+ * Returns a Voronoi diagram constructed
+ * from the points of the input geometry.
+ *
+ ******************************************/
+Datum ST_Voronoi(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_Voronoi);
+Datum ST_Voronoi(PG_FUNCTION_ARGS)
+{
+#if POSTGIS_GEOS_VERSION < 35
+ lwpgerror("The GEOS version this PostGIS binary "
+ "was compiled against (%d) doesn't support "
+ "'ST_Voronoi' function (3.5.0+ required)",
+ POSTGIS_GEOS_VERSION);
+ PG_RETURN_NULL();
+#else /* POSTGIS_GEOS_VERSION >= 35 */
+ GSERIALIZED* input;
+ GSERIALIZED* clip;
+ GSERIALIZED* result;
+ LWGEOM* lwgeom_input;
+ LWGEOM* lwgeom_result;
+ double tolerance;
+ GBOX clip_envelope;
+ int custom_clip_envelope;
+ int return_polygons;
+
+ /* Return NULL on NULL geometry */
+ if (PG_ARGISNULL(0))
+ PG_RETURN_NULL();
+
+ /* Read our tolerance value */
+ if (PG_ARGISNULL(2))
+ {
+ lwpgerror("Tolerance must be a positive number.");
+ PG_RETURN_NULL();
+ }
+
+ tolerance = PG_GETARG_FLOAT8(2);
+
+ if (tolerance < 0)
+ {
+ lwpgerror("Tolerance must be a positive number.");
+ PG_RETURN_NULL();
+ }
+
+ /* Are we returning lines or polygons? */
+ if (PG_ARGISNULL(3))
+ {
+ lwpgerror("return_polygons must be true or false.");
+ PG_RETURN_NULL();
+ }
+ return_polygons = PG_GETARG_BOOL(3);
+
+ /* Read our clipping envelope, if applicable. */
+ custom_clip_envelope = !PG_ARGISNULL(1);
+ if (custom_clip_envelope) {
+ clip = PG_GETARG_GSERIALIZED_P(1);
+ if (!gserialized_get_gbox_p(clip, &clip_envelope))
+ {
+ lwpgerror("Could not determine envelope of clipping geometry.");
+ PG_FREE_IF_COPY(clip, 0);
+ PG_RETURN_NULL();
+ }
+ PG_FREE_IF_COPY(clip, 0);
+ }
+
+ /* Read our input geometry */
+ input = PG_GETARG_GSERIALIZED_P(0);
+
+ lwgeom_input = lwgeom_from_gserialized(input);
+
+ if(!lwgeom_input)
+ {
+ lwpgerror("Could not read input geometry.");
+ PG_FREE_IF_COPY(input, 0);
+ PG_RETURN_NULL();
+ }
+
+ lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
+ lwgeom_free(lwgeom_input);
+
+ if (!lwgeom_result)
+ {
+ lwpgerror("Error computing Voronoi diagram.");
+ PG_FREE_IF_COPY(input, 0);
+ PG_RETURN_NULL();
+ }
+
+ result = geometry_serialize(lwgeom_result);
+ lwgeom_free(lwgeom_result);
+
+ PG_FREE_IF_COPY(input, 0);
+ PG_RETURN_POINTER(result);
+
+#endif /* POSTGIS_GEOS_VERSION >= 35 */
+}
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/postgis/postgis.sql.in 2015-12-02 01:10:33 UTC (rev 14467)
@@ -3516,7 +3516,37 @@
LANGUAGE 'c' IMMUTABLE STRICT
COST 100;
+--------------------------------------------------------------------------------
+-- ST_Voronoi
+--------------------------------------------------------------------------------
+-- ST_Voronoi(g1 geometry, clip geometry, tolerance float8, return_polygons boolean)
+--
+-- Builds a Voronoi Diagram from the vertices of the supplied geometry.
+--
+-- By default, the diagram will be extended to an envelope larger than the
+-- input points.
+--
+-- If a second geometry is supplied, the diagram will be extended to fill the
+-- envelope of the second geometry, unless that is smaller than the default
+-- envelope.
+--
+-- If a tolerance is given it will be used to snap the input points
+-- each-other.
+--
+-- If return_polygons is true, returns a GeometryCollection of polygons.
+-- If return_polygons is false, returns a MultiLineString.
+--
+-- Availability: 2.3.0
+-- Requires GEOS >= 3.5.0
+--
+CREATE OR REPLACE FUNCTION ST_Voronoi(g1 geometry, clip geometry DEFAULT NULL, tolerance float8 DEFAULT 0.0, return_polygons boolean DEFAULT true)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_Voronoi'
+ LANGUAGE 'c' IMMUTABLE
+ COST 100;
+
+
--------------------------------------------------------------------------------
-- Aggregates and their supporting functions
--------------------------------------------------------------------------------
Modified: trunk/regress/Makefile.in
===================================================================
--- trunk/regress/Makefile.in 2015-12-01 23:16:41 UTC (rev 14466)
+++ trunk/regress/Makefile.in 2015-12-02 01:10:33 UTC (rev 14467)
@@ -207,10 +207,11 @@
ifeq ($(shell expr $(POSTGIS_GEOS_VERSION) ">=" 35),1)
# GEOS-3.5 adds:
- # ST_ClipByBox2d
+ # ST_ClipByBox2d, ST_Subdivide, ST_Voronoi
TESTS += \
clipbybox2d \
- subdivide
+ subdivide \
+ voronoi
endif
Added: trunk/regress/voronoi.sql
===================================================================
--- trunk/regress/voronoi.sql (rev 0)
+++ trunk/regress/voronoi.sql 2015-12-02 01:10:33 UTC (rev 14467)
@@ -0,0 +1,17 @@
+-- postgres
+
+-- SRID is preserved
+SELECT 1, 32145 = ST_SRID(ST_Voronoi('SRID=32145;MULTIPOINT (0 0, 1 1, 2 2)'));
+-- NULL -> NULL
+SELECT 2, ST_Voronoi(NULL) IS NULL;
+-- NULL tolerance produces error
+SELECT 3, ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)', NULL, NULL);
+-- NULL return_polygons produces error
+SELECT 4, ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)', NULL, 0, NULL);
+-- Tolerance can't be negative
+SELECT 5, ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)', NULL, -2);
+-- Output types are correct
+SELECT 6, GeometryType(ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)')) = 'GEOMETRYCOLLECTION';
+SELECT 7, GeometryType(ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)', NULL, 0, false)) = 'MULTILINESTRING';
+-- Clipping extent is handled correctly
+SELECT 8, ST_Equals(ST_Envelope('LINESTRING (-20 -10, 10 10)'::geometry), ST_Envelope(ST_Voronoi('MULTIPOINT (0 0, 1 1, 2 2)', 'MULTIPOINT (-20 -10, 10 10)')));
Added: trunk/regress/voronoi_expected
===================================================================
--- trunk/regress/voronoi_expected (rev 0)
+++ trunk/regress/voronoi_expected 2015-12-02 01:10:33 UTC (rev 14467)
@@ -0,0 +1,8 @@
+1|t
+2|t
+ERROR: Tolerance must be a positive number.
+ERROR: return_polygons must be true or false.
+ERROR: Tolerance must be a positive number.
+6|t
+7|t
+8|t
More information about the postgis-tickets
mailing list