[postgis-tickets] r14449 - Add missing files for #2996
Daniel Baston
dbaston at gmail.com
Sun Nov 29 15:16:30 PST 2015
Author: dbaston
Date: 2015-11-29 15:16:30 -0800 (Sun, 29 Nov 2015)
New Revision: 14449
Added:
trunk/liblwgeom/cunit/cu_minimum_bounding_circle.c
trunk/liblwgeom/lwboundingcircle.c
trunk/regress/minimum_bounding_circle.sql
trunk/regress/minimum_bounding_circle_expected
Log:
Add missing files for #2996
Added: trunk/liblwgeom/cunit/cu_minimum_bounding_circle.c
===================================================================
--- trunk/liblwgeom/cunit/cu_minimum_bounding_circle.c (rev 0)
+++ trunk/liblwgeom/cunit/cu_minimum_bounding_circle.c 2015-11-29 23:16:30 UTC (rev 14449)
@@ -0,0 +1,80 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ * Copyright 2015 Daniel Baston
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "CUnit/Basic.h"
+#include "cu_tester.h"
+#include "../liblwgeom_internal.h"
+
+static void mbc_test(LWGEOM* g)
+{
+ LWBOUNDINGCIRCLE* result = lwgeom_calculate_mbc(g);
+ CU_ASSERT_TRUE(result != NULL);
+
+ LWPOINTITERATOR* it = lwpointiterator_create(g);
+
+ POINT2D p;
+ POINT4D p4;
+ while (lwpointiterator_next(it, &p4))
+ {
+ p.x = p4.x;
+ p.y = p4.y;
+
+ CU_ASSERT_TRUE(distance2d_pt_pt(result->center, &p) <= result->radius);
+ }
+
+ lwboundingcircle_destroy(result);
+ lwpointiterator_destroy(it);
+}
+
+static void basic_test(void)
+{
+ uint32_t i;
+
+ char* inputs[] =
+ {
+ "POLYGON((26426 65078,26531 65242,26075 65136,26096 65427,26426 65078))",
+ "POINT (17 253)",
+ "TRIANGLE ((0 0, 10 0, 10 10, 0 0))",
+ "LINESTRING (17 253, -44 28, 33 11, 26 44)",
+ "MULTIPOINT ((0 0), (0 0), (0 0), (0 0))",
+ "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))",
+ "LINESTRING (-48546889 37039202, -37039202 -48546889)"
+ };
+
+ for (i = 0; i < sizeof(inputs)/sizeof(LWGEOM*); i++)
+ {
+ LWGEOM* input = lwgeom_from_wkt(inputs[i], LW_PARSER_CHECK_NONE);
+ mbc_test(input);
+ lwgeom_free(input);
+ }
+
+}
+
+static void test_empty(void)
+{
+ LWGEOM* input = lwgeom_from_wkt("POINT EMPTY", LW_PARSER_CHECK_NONE);
+
+ LWBOUNDINGCIRCLE* result = lwgeom_calculate_mbc(input);
+ CU_ASSERT_TRUE(result == NULL);
+
+ lwgeom_free(input);
+}
+
+/*
+ ** Used by test harness to register the tests in this file.
+ */
+void minimum_bounding_circle_suite_setup(void);
+void minimum_bounding_circle_suite_setup(void)
+{
+ CU_pSuite suite = CU_add_suite("minimum_bounding_circle", NULL, NULL);
+ PG_ADD_TEST(suite, basic_test);
+ PG_ADD_TEST(suite, test_empty);
+}
Added: trunk/liblwgeom/lwboundingcircle.c
===================================================================
--- trunk/liblwgeom/lwboundingcircle.c (rev 0)
+++ trunk/liblwgeom/lwboundingcircle.c 2015-11-29 23:16:30 UTC (rev 14449)
@@ -0,0 +1,280 @@
+/**********************************************************************
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.net
+ *
+ * Copyright 2015 Daniel Baston <dbaston at gmail.com>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include <string.h>
+#include "liblwgeom_internal.h"
+
+typedef struct {
+ const POINT2D* p1;
+ const POINT2D* p2;
+ const POINT2D* p3;
+} SUPPORTING_POINTS;
+
+static SUPPORTING_POINTS*
+supporting_points_create()
+{
+ SUPPORTING_POINTS* s = lwalloc(sizeof(SUPPORTING_POINTS));
+ s->p1 = NULL;
+ s->p2 = NULL;
+ s->p3 = NULL;
+
+ return s;
+}
+
+static void
+supporting_points_destroy(SUPPORTING_POINTS* s)
+{
+ lwfree(s);
+}
+
+static uint32_t
+num_supporting_points(SUPPORTING_POINTS* support)
+{
+ uint32_t N = 0;
+
+ if (support->p1 != NULL)
+ N++;
+ if (support->p2 != NULL)
+ N++;
+ if (support->p3 != NULL)
+ N++;
+
+ return N;
+}
+
+static int
+add_supporting_point(SUPPORTING_POINTS* support, const POINT2D* p)
+{
+ switch(num_supporting_points(support))
+ {
+ case 0: support->p1 = p;
+ break;
+ case 1: support->p2 = p;
+ break;
+ case 2: support->p3 = p;
+ break;
+ default: return LW_FAILURE;
+ }
+
+ return LW_SUCCESS;
+}
+
+static int
+point_inside_circle(const POINT2D* p, const LWBOUNDINGCIRCLE* c)
+{
+ if (!c)
+ return LW_FALSE;
+
+ if (distance2d_pt_pt(p, c->center) > c->radius)
+ return LW_FALSE;
+
+ return LW_TRUE;
+}
+
+static double
+det(double m00, double m01, double m10, double m11)
+{
+ return m00 * m11 - m01 * m10;
+}
+
+static void
+circumcenter(const POINT2D* a, const POINT2D* b, const POINT2D* c, POINT2D* result)
+{
+ double cx = c->x;
+ double cy = c->y;
+ double ax = a->x - cx;
+ double ay = a->y - cy;
+ double bx = b->x - cx;
+ double by = b->y - cy;
+
+ double denom = 2 * det(ax, ay, bx, by);
+ double numx = det(ay, ax * ax + ay * ay, by, bx * bx + by * by);
+ double numy = det(ax, ax * ax + ay * ay, bx, bx * bx + by * by);
+
+ result->x = cx - numx / denom;
+ result->y = cy + numy / denom;
+}
+
+static void
+calculate_mbc_1(const SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
+{
+ mbc->radius = 0;
+ mbc->center->x = support->p1->x;
+ mbc->center->y = support->p1->y;
+}
+
+static void
+calculate_mbc_2(const SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
+{
+ double d1, d2;
+
+ mbc->center->x = 0.5*(support->p1->x + support->p2->x);
+ mbc->center->y = 0.5*(support->p1->y + support->p2->y);
+
+ d1 = distance2d_pt_pt(mbc->center, support->p1);
+ d2 = distance2d_pt_pt(mbc->center, support->p2);
+
+ mbc->radius = FP_MAX(d1, d2);
+}
+
+static void
+calculate_mbc_3(const SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
+{
+ double d1, d2, d3;
+ circumcenter(support->p1, support->p2, support->p3, mbc->center);
+
+ d1 = distance2d_pt_pt(mbc->center, support->p1);
+ d2 = distance2d_pt_pt(mbc->center, support->p2);
+ d3 = distance2d_pt_pt(mbc->center, support->p3);
+
+ mbc->radius = FP_MAX(FP_MAX(d1, d2), d3);
+}
+
+static int
+calculate_mbc_from_support(SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
+{
+ switch(num_supporting_points(support))
+ {
+ case 0: break;
+ case 1: calculate_mbc_1(support, mbc);
+ break;
+ case 2: calculate_mbc_2(support, mbc);
+ break;
+ case 3: calculate_mbc_3(support, mbc);
+ break;
+ default: return LW_FAILURE;
+ }
+
+ return LW_SUCCESS;
+}
+
+static int
+calculate_mbc(const POINT2D** points, uint32_t max_n, SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
+{
+ uint32_t i;
+
+ if(!calculate_mbc_from_support(support, mbc))
+ {
+ return LW_FAILURE;
+ }
+
+ if (num_supporting_points(support) == 3)
+ {
+ /* If we're entering the function with three supporting points already, our circle
+ * is already fully constrained - we couldn't add another supporting point if we
+ * needed to. So, there's no point in proceeding further. Welzl (1991) provides
+ * a much better explanation of why this works.
+ * */
+ return LW_SUCCESS;
+ }
+
+ for (i = 0; i < max_n; i++)
+ {
+ if (!point_inside_circle(points[i], mbc))
+ {
+ /* We've run into a point that isn't inside our circle. To fix this, we'll
+ * go back in time, and re-run the algorithm for each point we've seen so
+ * far, with the constraint that the current point must be on the boundary
+ * of the circle. Then, we'll continue on in this loop with the modified
+ * circle that by definition includes the current point. */
+ SUPPORTING_POINTS next_support;
+ memcpy(&next_support, support, sizeof(SUPPORTING_POINTS));
+
+ add_supporting_point(&next_support, points[i]);
+ if (!calculate_mbc(points, i, &next_support, mbc))
+ {
+ return LW_FAILURE;
+ }
+ }
+ }
+
+ return LW_SUCCESS;
+}
+
+static LWBOUNDINGCIRCLE*
+lwboundingcircle_create()
+{
+ LWBOUNDINGCIRCLE* c = lwalloc(sizeof(LWBOUNDINGCIRCLE));
+ c->center = lwalloc(sizeof(POINT2D));
+
+ c->radius = 0.0;
+ c->center->x = 0.0;
+ c->center->y = 0.0;
+
+ return c;
+}
+
+void
+lwboundingcircle_destroy(LWBOUNDINGCIRCLE* c)
+{
+ lwfree(c->center);
+ lwfree(c);
+}
+
+LWBOUNDINGCIRCLE*
+lwgeom_calculate_mbc(const LWGEOM* g)
+{
+ SUPPORTING_POINTS* support;
+ LWBOUNDINGCIRCLE* result;
+ LWPOINTITERATOR* it;
+ uint32_t num_points;
+ POINT2D** points;
+ POINT4D p;
+ uint32_t i;
+ int success;
+
+ if(g == NULL || lwgeom_is_empty(g))
+ return LW_FAILURE;
+
+ num_points = lwgeom_count_vertices(g);
+ it = lwpointiterator_create(g);
+ points = lwalloc(num_points * sizeof(POINT2D*));
+ for (i = 0; i < num_points; i++)
+ {
+ if(!lwpointiterator_next(it, &p))
+ {
+ uint32_t j;
+ for (j = 0; j < i; j++)
+ {
+ lwfree(points[j]);
+ }
+ lwpointiterator_destroy(it);
+ lwfree(points);
+ return LW_FAILURE;
+ }
+
+ points[i] = lwalloc(sizeof(POINT2D));
+ points[i]->x = p.x;
+ points[i]->y = p.y;
+ }
+ lwpointiterator_destroy(it);
+
+ support = supporting_points_create();
+ result = lwboundingcircle_create();
+ /* Technically, a randomized algorithm would demand that we shuffle the input points
+ * before we call calculate_mbc(). However, we make the (perhaps poor) assumption that
+ * the order we happen to find the points is as good as random, or close enough.
+ * */
+ success = calculate_mbc((const POINT2D**) points, num_points, support, result);
+
+ for (i = 0; i < num_points; i++)
+ {
+ lwfree(points[i]);
+ }
+ lwfree(points);
+ supporting_points_destroy(support);
+
+ if (!success)
+ return NULL;
+
+ return result;
+}
Added: trunk/regress/minimum_bounding_circle.sql
===================================================================
--- trunk/regress/minimum_bounding_circle.sql (rev 0)
+++ trunk/regress/minimum_bounding_circle.sql 2015-11-29 23:16:30 UTC (rev 14449)
@@ -0,0 +1,6 @@
+SELECT 't1', ST_MinimumBoundingCircle(NULL::geometry) IS NULL;
+SELECT 't2', ST_Equals(ST_MinimumBoundingCircle('POINT EMPTY'), 'POLYGON EMPTY'::geometry);
+SELECT 't3', ST_SRID(ST_MinimumBoundingCircle('SRID=32611;POINT(4021690.58034526 6040138.01373556)')) = 32611;
+SELECT 't4', ST_SRID(center) = 32611 FROM ST_MinimumBoundingRadius('SRID=32611;POINT(4021690.58034526 6040138.01373556)');
+SELECT 't5', ST_Equals(center, 'POINT EMPTY') AND radius = 0 FROM ST_MinimumBoundingRadius('GEOMETRYCOLLECTION EMPTY');
+SELECT 't6', ST_Equals(center, 'POINT (0 0.5)') AND radius = 0.5 FROM ST_MinimumBoundingRadius('MULTIPOINT ((0 0 0), (0 1 1))');
Added: trunk/regress/minimum_bounding_circle_expected
===================================================================
--- trunk/regress/minimum_bounding_circle_expected (rev 0)
+++ trunk/regress/minimum_bounding_circle_expected 2015-11-29 23:16:30 UTC (rev 14449)
@@ -0,0 +1,6 @@
+t1|t
+t2|t
+t3|t
+t4|t
+t5|t
+t6|t
More information about the postgis-tickets
mailing list