New method to approximate minimum bounding circle polygon without using ST_Buffer...

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 da44104c145b42fd9ae9330d8bb43e5f6ebd8846..aed358fc6fbdb145e6eba25a4105b950aed16ef4 100644 (file)
--- 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 adb34dbb303fb5631ed1583d5484e94b81b52a4d..e96839e10cc8a537d339ce16208efbf42d50b7cf 100644 (file)
--- 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 115c0a5dd576ed3e8cc9c5b46a6c4be835102ad6..293133060d5e472dac7f1367b7166dd555c80d76 100644 (file)
--- a/liblwgeom/lwpoly.c
+++ b/liblwgeom/lwpoly.c
@@ -29,6 +29,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 #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 c1dbb97cf57054de85aa13308673ae5841845f84..922f5d98b103bd4808d4b538ed772c608330bf0f 100644 (file)
--- 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 d97f82bb907d5da89c04d69a603cab305d569bb4..8c33224d4ac92be8e2787d1e4f3b5466df071601 100644 (file)
--- 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 4c6cb5ce28828e2d9b6609621f43e17838c238fe..9672c6759f20c9496030e7080ab4b0ccb66e3a6b 100644 (file)
--- 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 c2a59c2d9882f6c7287746edb13d89510b864a66..b5bdf4814ddbb0d7b814f90333e226e395bde5d4 100644 (file)
--- 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