From: Daniel Baston Date: Wed, 14 Sep 2016 16:04:46 +0000 (+0000) Subject: New method to approximate minimum bounding circle polygon without using ST_Buffer... X-Git-Tag: 2.3.0rc1~7 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=bd3177743db5cb786deb05d7327f39a5cec2daf6;p=postgis New method to approximate minimum bounding circle polygon without using ST_Buffer. (References #2841, #3620) git-svn-id: http://svn.osgeo.org/postgis/trunk@15109 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c index da44104c1..aed358fc6 100644 --- a/liblwgeom/cunit/cu_algorithm.c +++ b/liblwgeom/cunit/cu_algorithm.c @@ -1102,6 +1102,46 @@ static void test_point_density(void) 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 @@ void algorithms_suite_setup(void) 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); } diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in index adb34dbb3..e96839e10 100644 --- a/liblwgeom/liblwgeom.h.in +++ b/liblwgeom/liblwgeom.h.in @@ -1357,6 +1357,7 @@ extern LWLINE *lwline_removepoint(LWLINE *line, uint32_t which); 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 */ diff --git a/liblwgeom/lwpoly.c b/liblwgeom/lwpoly.c index 115c0a5dd..293133060 100644 --- a/liblwgeom/lwpoly.c +++ b/liblwgeom/lwpoly.c @@ -29,6 +29,7 @@ #include #include #include +#include #include "liblwgeom_internal.h" #include "lwgeom_log.h" @@ -93,6 +94,45 @@ lwpoly_construct_rectangle(char hasz, char hasm, POINT4D *p1, POINT4D *p2, return lwpoly; } +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) { diff --git a/postgis/lwgeom_functions_analytic.c b/postgis/lwgeom_functions_analytic.c index c1dbb97cf..922f5d98b 100644 --- a/postgis/lwgeom_functions_analytic.c +++ b/postgis/lwgeom_functions_analytic.c @@ -53,6 +53,7 @@ Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS); 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 @@ Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS) 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); @@ -1142,6 +1143,61 @@ Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS) PG_RETURN_DATUM(result); } +/********************************************************************** + * + * 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 diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in index d97f82bb9..8c33224d4 100644 --- a/postgis/postgis.sql.in +++ b/postgis/postgis.sql.in @@ -3305,9 +3305,9 @@ CREATE OR REPLACE FUNCTION ST_MinimumBoundingRadius(geometry, OUT center geometr -- Availability: 1.4.0 CREATE OR REPLACE FUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer DEFAULT 48) - RETURNS geometry - AS $$ SELECT @extschema@.ST_Buffer(center, radius, segs_per_quarter) FROM @extschema@.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 '') diff --git a/regress/tickets.sql b/regress/tickets.sql index 4c6cb5ce2..9672c6759 100644 --- a/regress/tickets.sql +++ b/regress/tickets.sql @@ -997,6 +997,10 @@ select 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); diff --git a/regress/tickets_expected b/regress/tickets_expected index c2a59c2d9..b5bdf4814 100644 --- a/regress/tickets_expected +++ b/regress/tickets_expected @@ -296,5 +296,7 @@ ERROR: invalid KML representation #3470b|50 #3569|BOX(1 2,3 4) #3579|f|t +#3620a|POINT(3 7) +#3620b|POINT(2 8) #3627a|o}~~|AdshNoSsBgd@eGoBlm@wKhj@~@? #3627b|t