[postgis-tickets] r15109 - New method to approximate minimum bounding circle polygon without using ST_Buffer. (References #2841, #3620)
Daniel Baston
dbaston at gmail.com
Wed Sep 14 09:04:47 PDT 2016
Author: dbaston
Date: 2016-09-14 09:04:46 -0700 (Wed, 14 Sep 2016)
New Revision: 15109
Modified:
trunk/liblwgeom/cunit/cu_algorithm.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/lwpoly.c
trunk/postgis/lwgeom_functions_analytic.c
trunk/postgis/postgis.sql.in
trunk/regress/tickets.sql
trunk/regress/tickets_expected
Log:
New method to approximate minimum bounding circle polygon without using ST_Buffer. (References #2841, #3620)
Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/liblwgeom/cunit/cu_algorithm.c 2016-09-14 16:04:46 UTC (rev 15109)
@@ -1102,6 +1102,46 @@
lwgeom_free(geom);
}
+static void test_lwpoly_construct_circle(void)
+{
+ LWPOLY* p;
+ GBOX* g;
+ const int srid = 4326;
+ const int segments_per_quad = 17;
+ const int x = 10;
+ const int y = 20;
+ const int r = 5;
+
+ /* With normal arguments you should get something circle-y */
+ p = lwpoly_construct_circle(srid, x, y, r, segments_per_quad, LW_TRUE);
+
+ ASSERT_INT_EQUAL(lwgeom_count_vertices(p), segments_per_quad * 4 + 1);
+ ASSERT_INT_EQUAL(lwgeom_get_srid(lwpoly_as_lwgeom(p)), srid)
+
+ g = lwgeom_get_bbox(lwpoly_as_lwgeom(p));
+ CU_ASSERT_DOUBLE_EQUAL(g->xmin, x-r, 0.1);
+ CU_ASSERT_DOUBLE_EQUAL(g->xmax, x+r, 0.1);
+ CU_ASSERT_DOUBLE_EQUAL(g->ymin, y-r, 0.1);
+ CU_ASSERT_DOUBLE_EQUAL(g->ymax, y+r, 0.1);
+
+ CU_ASSERT_DOUBLE_EQUAL(lwgeom_area(p), M_PI*5*5, 0.1);
+
+ lwpoly_free(p);
+
+ /* No segments? No circle. */
+ p = lwpoly_construct_circle(srid, x, y, r, 0, LW_TRUE);
+ CU_ASSERT_TRUE(p == NULL);
+
+ /* Negative radius? No circle. */
+ p = lwpoly_construct_circle(srid, x, y, -1e-3, segments_per_quad, LW_TRUE);
+ CU_ASSERT_TRUE(p == NULL);
+
+ /* Zero radius? Invalid circle */
+ p = lwpoly_construct_circle(srid, x, y, 0, segments_per_quad, LW_TRUE);
+ CU_ASSERT_TRUE(p != NULL);
+ lwpoly_free(p);
+}
+
static void test_kmeans(void)
{
static int cluster_size = 100;
@@ -1172,4 +1212,5 @@
PG_ADD_TEST(suite,test_kmeans);
PG_ADD_TEST(suite,test_median_handles_3d_correctly);
PG_ADD_TEST(suite,test_median_robustness);
+ PG_ADD_TEST(suite,test_lwpoly_construct_circle);
}
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/liblwgeom/liblwgeom.h.in 2016-09-14 16:04:46 UTC (rev 15109)
@@ -1357,6 +1357,7 @@
extern void lwline_setPoint4d(LWLINE *line, uint32_t which, POINT4D *newpoint);
extern LWPOLY *lwpoly_from_lwlines(const LWLINE *shell, uint32_t nholes, const LWLINE **holes);
extern LWPOLY* lwpoly_construct_rectangle(char hasz, char hasm, POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D *p4);
+extern LWPOLY* lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior);
extern LWTRIANGLE *lwtriangle_from_lwline(const LWLINE *shell);
extern LWMPOINT *lwmpoint_from_lwgeom(const LWGEOM *g); /* Extract the coordinates of an LWGEOM into an LWMPOINT */
Modified: trunk/liblwgeom/lwpoly.c
===================================================================
--- trunk/liblwgeom/lwpoly.c 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/liblwgeom/lwpoly.c 2016-09-14 16:04:46 UTC (rev 15109)
@@ -29,6 +29,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
@@ -94,6 +95,45 @@
}
LWPOLY*
+lwpoly_construct_circle(int srid, double x, double y, double radius, uint32_t segments_per_quarter, char exterior)
+{
+ const int segments = 4*segments_per_quarter;
+ const double theta = 2*M_PI / segments;
+ LWPOLY *lwpoly;
+ POINTARRAY *pa;
+ POINT4D pt;
+ uint32_t i;
+
+ if (segments_per_quarter < 1)
+ {
+ lwerror("Need at least one segment per quarter-circle.");
+ return NULL;
+ }
+
+ if (radius < 0)
+ {
+ lwerror("Radius must be positive.");
+ return NULL;
+ }
+
+ lwpoly = lwpoly_construct_empty(srid, LW_FALSE, LW_FALSE);
+ pa = ptarray_construct_empty(LW_FALSE, LW_FALSE, segments + 1);
+
+ if (exterior)
+ radius *= sqrt(1 + pow(tan(theta/2), 2));
+
+ for (i = 0; i <= segments; i++)
+ {
+ pt.x = x + radius*sin(i * theta);
+ pt.y = y + radius*cos(i * theta);
+ ptarray_append_point(pa, &pt, LW_TRUE);
+ }
+
+ lwpoly_add_ring(lwpoly, pa);
+ return lwpoly;
+}
+
+LWPOLY*
lwpoly_construct_empty(int srid, char hasz, char hasm)
{
LWPOLY *result = lwalloc(sizeof(LWPOLY));
Modified: trunk/postgis/lwgeom_functions_analytic.c
===================================================================
--- trunk/postgis/lwgeom_functions_analytic.c 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/postgis/lwgeom_functions_analytic.c 2016-09-14 16:04:46 UTC (rev 15109)
@@ -53,6 +53,7 @@
Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS);
Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS);
+Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS);
Datum ST_GeometricMedian(PG_FUNCTION_ARGS);
@@ -1110,7 +1111,7 @@
input = lwgeom_from_gserialized(geom);
mbc = lwgeom_calculate_mbc(input);
- if (!mbc)
+ if (!(mbc && mbc->center))
{
lwpgerror("Error calculating minimum bounding circle.");
lwgeom_free(input);
@@ -1144,6 +1145,61 @@
/**********************************************************************
*
+ * ST_MinimumBoundingCircle
+ *
+ **********************************************************************/
+
+PG_FUNCTION_INFO_V1(ST_MinimumBoundingCircle);
+Datum ST_MinimumBoundingCircle(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED* geom;
+ LWGEOM* input;
+ LWBOUNDINGCIRCLE* mbc = NULL;
+ LWGEOM* lwcircle;
+ GSERIALIZED* center;
+ int segs_per_quarter;
+
+ if (PG_ARGISNULL(0))
+ PG_RETURN_NULL();
+
+ geom = PG_GETARG_GSERIALIZED_P(0);
+ segs_per_quarter = PG_GETARG_INT32(1);
+
+ /* Empty geometry? Return POINT EMPTY */
+ if (gserialized_is_empty(geom))
+ {
+ lwcircle = (LWGEOM*) lwpoint_construct_empty(gserialized_get_srid(geom), LW_FALSE, LW_FALSE);
+ }
+ else
+ {
+ input = lwgeom_from_gserialized(geom);
+ mbc = lwgeom_calculate_mbc(input);
+
+ if (!(mbc && mbc->center))
+ {
+ lwpgerror("Error calculating minimum bounding circle.");
+ lwgeom_free(input);
+ PG_RETURN_NULL();
+ }
+
+ /* Zero radius? Return a point. */
+ if (mbc->radius == 0)
+ lwcircle = lwpoint_as_lwgeom(lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y));
+ else
+ lwcircle = lwpoly_as_lwgeom(lwpoly_construct_circle(input->srid, mbc->center->x, mbc->center->y, mbc->radius, segs_per_quarter, LW_TRUE));
+
+ lwboundingcircle_destroy(mbc);
+ lwgeom_free(input);
+ }
+
+ center = geometry_serialize(lwcircle);
+ lwgeom_free(lwcircle);
+
+ PG_RETURN_POINTER(center);
+}
+
+/**********************************************************************
+ *
* ST_GeometricMedian
*
**********************************************************************/
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/postgis/postgis.sql.in 2016-09-14 16:04:46 UTC (rev 15109)
@@ -3305,9 +3305,9 @@
-- Availability: 1.4.0
CREATE OR REPLACE FUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer DEFAULT 48)
- RETURNS geometry
- AS $$ SELECT @extschema at .ST_Buffer(center, radius, segs_per_quarter) FROM @extschema at .ST_MinimumBoundingRadius($1) sq $$
- LANGUAGE 'sql' IMMUTABLE STRICT _PARALLEL;
+ RETURNS geometry
+ AS 'MODULE_PATHNAME', 'ST_MinimumBoundingCircle'
+ LANGUAGE 'c' IMMUTABLE STRICT _PARALLEL;
-- Availability: 2.0.0 - requires GEOS-3.2 or higher
CREATE OR REPLACE FUNCTION ST_OffsetCurve(line geometry, distance float8, params text DEFAULT '')
Modified: trunk/regress/tickets.sql
===================================================================
--- trunk/regress/tickets.sql 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/regress/tickets.sql 2016-09-14 16:04:46 UTC (rev 15109)
@@ -997,6 +997,10 @@
ST_Intersects(ST_Buffer(road.geom, sidewalk_offset + epsilon), sidewalks.geom) -- should be true
from road, sidewalks, params;
+-- #3620
+SELECT '#3620a', ST_AsText(ST_MinimumBoundingCircle('POINT (3 7)'));
+SELECT '#3620b', ST_AsText(ST_MinimumBoundingCircle('LINESTRING (2 8, 2 8)'));
+
-- #3627
SELECT '#3627a', ST_AsEncodedPolyline('SRID=4326;LINESTRING(-0.250691 49.283048,-0.250633 49.283376,-0.250502 49.283972,-0.251245 49.284028,-0.251938 49.284232,-0.251938 49.2842)', 6);
SELECT '#3627b', ST_Equals(geom, ST_LineFromEncodedPolyline(ST_AsEncodedPolyline(geom, 7), 7)) FROM (VALUES ('SRID=4326;LINESTRING (0 0, 1 1)')) AS v (geom);
Modified: trunk/regress/tickets_expected
===================================================================
--- trunk/regress/tickets_expected 2016-09-14 13:37:28 UTC (rev 15108)
+++ trunk/regress/tickets_expected 2016-09-14 16:04:46 UTC (rev 15109)
@@ -296,5 +296,7 @@
#3470b|50
#3569|BOX(1 2,3 4)
#3579|f|t
+#3620a|POINT(3 7)
+#3620b|POINT(2 8)
#3627a|o}~~|AdshNoSsBgd at eGoBlm@wKhj@~@?
#3627b|t
More information about the postgis-tickets
mailing list