[postgis-tickets] r14448 - #2996 ST_MinimumBoundingCircle doesn't always contain all points
Daniel Baston
dbaston at gmail.com
Sun Nov 29 15:10:49 PST 2015
Author: dbaston
Date: 2015-11-29 15:10:49 -0800 (Sun, 29 Nov 2015)
New Revision: 14448
Modified:
trunk/doc/reference_processing.xml
trunk/liblwgeom/Makefile.in
trunk/liblwgeom/cunit/Makefile.in
trunk/liblwgeom/cunit/cu_tester.c
trunk/liblwgeom/liblwgeom.h.in
trunk/postgis/lwgeom_functions_analytic.c
trunk/postgis/lwgeom_geos.c
trunk/postgis/postgis.sql.in
trunk/regress/Makefile.in
trunk/regress/tickets.sql
trunk/regress/tickets_expected
Log:
#2996 ST_MinimumBoundingCircle doesn't always contain all points
Modified: trunk/doc/reference_processing.xml
===================================================================
--- trunk/doc/reference_processing.xml 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/doc/reference_processing.xml 2015-11-29 23:10:49 UTC (rev 14448)
@@ -2000,7 +2000,7 @@
<refsection>
<title>Description</title>
<para>Returns the smallest circle polygon that can fully contain a geometry. </para>
- <note><para>The circle is approximated by a polygon with a default of 48 segments per quarter circle. This number can be increased with little performance penalty to obtain a more accurate result.</para></note>
+ <note><para>The circle is approximated by a polygon with a default of 48 segments per quarter circle. Because the polygon is an approximation of the minimum bounding circle, some points in the input geometry may not be contained within the polygon. The approximation can be improved by increasing the number of segments, with little performance penalty. For applications where a polygonal approximation is not suitable, ST_MinimumBoundingRadius may be used.</para></note>
<para>It is often used with MULTI and Geometry Collections.
Although it is not an aggregate - you can use it in conjunction
@@ -2011,8 +2011,11 @@
<para>Availability: 1.4.0 - requires GEOS</para>
-
</refsection>
+ <refsection>
+ <title>See Also</title>
+ <para><xref linkend="ST_Collect" />, <xref linkend="ST_MinimumBoundingRadius" /></para>
+ </refsection>
<refsection>
<title>Examples</title>
@@ -2045,10 +2048,49 @@
</refsection>
<refsection>
<title>See Also</title>
- <para><xref linkend="ST_Collect" />, <xref linkend="ST_ConvexHull" /></para>
+ <para><xref linkend="ST_Collect" />, <xref linkend="ST_MinimumBoundingRadius" /></para>
</refsection>
</refentry>
+ <refentry id="ST_MinimumBoundingRadius">
+ <refnamediv>
+ <refname>ST_MinimumBoundingRadius</refname>
+ <refpurpose>Returns the center point and radius of the smallest circle that can fully contain a geometry.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>(geometry, double precision) <function>ST_MinimumBoundingRadius</function></funcdef>
+ <paramdef><type>geometry</type> <parameter>geom</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+ <para>Returns a record containing the center point and radius of the smallest circle that can fully contain a geometry.</para>
+ <para>Can be used in conjunction with <xref linkend="ST_Collect"/> to get the minimum bounding circle of a set of geometries.</para>
+ <para>Availability - 2.3.0</para>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+ <para><xref linkend="ST_Collect" />, <xref linkend="ST_MinimumBoundingCircle" /></para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+<programlisting>SELECT ST_AsText(center), radius FROM ST_MinimumBoundingRadius('POLYGON((26426 65078,26531 65242,26075 65136,26096 65427,26426 65078))');
+
+ st_astext | radius
+------------------------------------------+------------------
+ POINT(26284.8418027133 65267.1145090825) | 247.436045591407
+</programlisting>
+ </refsection>
+
+ </refentry>
+
<refentry id="ST_Polygonize">
<refnamediv>
<refname>ST_Polygonize</refname>
Modified: trunk/liblwgeom/Makefile.in
===================================================================
--- trunk/liblwgeom/Makefile.in 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/liblwgeom/Makefile.in 2015-11-29 23:10:49 UTC (rev 14448)
@@ -52,6 +52,7 @@
lwmpoint.o \
lwmline.o \
lwmpoly.o \
+ lwboundingcircle.o \
lwcollection.o \
lwcircstring.o \
lwcompound.o \
Modified: trunk/liblwgeom/cunit/Makefile.in
===================================================================
--- trunk/liblwgeom/cunit/Makefile.in 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/liblwgeom/cunit/Makefile.in 2015-11-29 23:10:49 UTC (rev 14448)
@@ -27,6 +27,7 @@
cu_buildarea.o \
cu_clean.o \
cu_print.o \
+ cu_minimum_bounding_circle.o \
cu_misc.o \
cu_ptarray.o \
cu_geodetic.o \
Modified: trunk/liblwgeom/cunit/cu_tester.c
===================================================================
--- trunk/liblwgeom/cunit/cu_tester.c 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/liblwgeom/cunit/cu_tester.c 2015-11-29 23:10:49 UTC (rev 14448)
@@ -43,6 +43,7 @@
extern void libgeom_suite_setup(void);
extern void measures_suite_setup(void);
extern void effectivearea_suite_setup(void);
+extern void minimum_bounding_circle_suite_setup(void);
extern void misc_suite_setup(void);
extern void node_suite_setup(void);
extern void out_encoded_polyline_suite_setup(void);
@@ -88,6 +89,7 @@
libgeom_suite_setup,
measures_suite_setup,
effectivearea_suite_setup,
+ minimum_bounding_circle_suite_setup,
misc_suite_setup,
node_suite_setup,
out_encoded_polyline_suite_setup,
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/liblwgeom/liblwgeom.h.in 2015-11-29 23:10:49 UTC (rev 14448)
@@ -1568,6 +1568,27 @@
*/
extern int lwgeom_covers_lwgeom_sphere(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2);
+typedef struct {
+ POINT2D* center;
+ double radius;
+} LWBOUNDINGCIRCLE;
+
+extern void lwboundingcircle_destroy(LWBOUNDINGCIRCLE* c);
+
+/* Calculates the minimum circle that encloses all of the points in g, using a
+ * two-dimensional implementation of the algorithm proposed in:
+ *
+ * Welzl, Emo (1991), "Smallest enclosing disks (balls and elipsoids)."
+ * New Results and Trends in Computer Science (H. Maurer, Ed.), Lecture Notes
+ * in Computer Science, 555 (1991) 359-370.
+ *
+ * Available online at the time of this writing at
+ * https://www.inf.ethz.ch/personal/emo/PublFiles/SmallEnclDisk_LNCS555_91.pdf
+ *
+ * Returns NULL if the circle could not be calculated.
+ */
+extern LWBOUNDINGCIRCLE* lwgeom_calculate_mbc(const LWGEOM* g);
+
/**
* Remove repeated points!
*/
Modified: trunk/postgis/lwgeom_functions_analytic.c
===================================================================
--- trunk/postgis/lwgeom_functions_analytic.c 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/postgis/lwgeom_functions_analytic.c 2015-11-29 23:10:49 UTC (rev 14448)
@@ -11,6 +11,7 @@
**********************************************************************/
#include "postgres.h"
+#include "funcapi.h"
#include "fmgr.h"
#include "liblwgeom.h"
#include "liblwgeom_internal.h" /* For FP comparators. */
@@ -19,6 +20,11 @@
#include "lwgeom_rtree.h"
#include "lwgeom_functions_analytic.h"
+#if POSTGIS_PGSQL_VERSION >= 93
+#include "access/htup_details.h"
+#else
+#include "access/htup.h"
+#endif
/***********************************************************************
* Simple Douglas-Peucker line simplification.
@@ -33,6 +39,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);
double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
@@ -1055,3 +1062,71 @@
* End of "Fast Winding Number Inclusion of a Point in a Polygon" derivative.
******************************************************************************/
+/**********************************************************************
+ *
+ * ST_MinimumBoundingRadius
+ *
+ **********************************************************************/
+
+PG_FUNCTION_INFO_V1(ST_MinimumBoundingRadius);
+Datum ST_MinimumBoundingRadius(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED* geom;
+ LWGEOM* input;
+ LWBOUNDINGCIRCLE* mbc = NULL;
+ LWGEOM* lwcenter;
+ GSERIALIZED* center;
+ TupleDesc resultTupleDesc;
+ HeapTuple resultTuple;
+ Datum result;
+ Datum result_values[2];
+ bool result_is_null[2];
+ double radius = 0;
+
+ if (PG_ARGISNULL(0))
+ PG_RETURN_NULL();
+
+ geom = PG_GETARG_GSERIALIZED_P(0);
+
+ /* Empty geometry? Return POINT EMPTY with zero radius */
+ if (gserialized_is_empty(geom))
+ {
+ lwcenter = (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)
+ {
+ lwpgerror("Error calculating minimum bounding circle.");
+ lwgeom_free(input);
+ PG_RETURN_NULL();
+ }
+
+ lwcenter = (LWGEOM*) lwpoint_make2d(input->srid, mbc->center->x, mbc->center->y);
+ radius = mbc->radius;
+
+ lwboundingcircle_destroy(mbc);
+ lwgeom_free(input);
+ }
+
+ center = geometry_serialize(lwcenter);
+ lwgeom_free(lwcenter);
+
+ get_call_result_type(fcinfo, NULL, &resultTupleDesc);
+ BlessTupleDesc(resultTupleDesc);
+
+ result_values[0] = PointerGetDatum(center);
+ result_is_null[0] = false;
+ result_values[1] = Float8GetDatum(radius);
+ result_is_null[1] = false;
+
+ resultTuple = heap_form_tuple(resultTupleDesc, result_values, result_is_null);
+
+ result = HeapTupleGetDatum(resultTuple);
+
+ PG_RETURN_DATUM(result);
+}
+
Modified: trunk/postgis/lwgeom_geos.c
===================================================================
--- trunk/postgis/lwgeom_geos.c 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/postgis/lwgeom_geos.c 2015-11-29 23:10:49 UTC (rev 14448)
@@ -3875,3 +3875,4 @@
#endif /* POSTGIS_GEOS_VERSION >= 33 */
}
+
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/postgis/postgis.sql.in 2015-11-29 23:10:49 UTC (rev 14448)
@@ -3174,6 +3174,17 @@
$$
LANGUAGE 'sql' IMMUTABLE STRICT;
+-- Availability: 2.3.0
+CREATE OR REPLACE FUNCTION ST_MinimumBoundingRadius(geometry, OUT center geometry, OUT radius double precision)
+ AS 'MODULE_PATHNAME', 'ST_MinimumBoundingRadius'
+ LANGUAGE 'c' IMMUTABLE STRICT;
+
+-- Availability: 1.4.0
+CREATE OR REPLACE FUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer DEFAULT 48)
+ RETURNS geometry
+ AS $$ SELECT ST_Buffer(center, radius, segs_per_quarter) FROM ST_MinimumBoundingRadius($1) sq $$
+ LANGUAGE 'sql' IMMUTABLE STRICT;
+
-- Availability: 2.0.0 - requires GEOS-3.2 or higher
CREATE OR REPLACE FUNCTION ST_OffsetCurve(line geometry, distance float8, params text DEFAULT '')
RETURNS geometry
@@ -5321,117 +5332,6 @@
-- USER CONTRIBUTED
---------------------------------------------------------------
------------------------------------------------------------------------
--- ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer)
------------------------------------------------------------------------
--- Returns the smallest circle polygon that can fully contain a geometry
--- Defaults to 48 segs per quarter to approximate a circle
--- Contributed by Bruce Rindahl
--- Availability: 1.4.0
------------------------------------------------------------------------
-CREATE OR REPLACE FUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer DEFAULT 48)
- RETURNS geometry AS
-$BODY$
- DECLARE
- hull GEOMETRY;
- ring GEOMETRY;
- center GEOMETRY;
- radius DOUBLE PRECISION;
- dist DOUBLE PRECISION;
- d DOUBLE PRECISION;
- idx1 integer;
- idx2 integer;
- l1 GEOMETRY;
- l2 GEOMETRY;
- p1 GEOMETRY;
- p2 GEOMETRY;
- a1 DOUBLE PRECISION;
- a2 DOUBLE PRECISION;
-
-
- BEGIN
-
- -- First compute the ConvexHull of the geometry
- hull = ST_ConvexHull(inputgeom);
- --A point really has no MBC
- IF ST_GeometryType(hull) = 'ST_Point' THEN
- RETURN hull;
- END IF;
- -- convert the hull perimeter to a linestring so we can manipulate individual points
- --If its already a linestring force it to a closed linestring
- ring = CASE WHEN ST_GeometryType(hull) = 'ST_LineString' THEN ST_AddPoint(hull, ST_StartPoint(hull)) ELSE ST_ExteriorRing(hull) END;
-
- dist = 0;
- -- Brute Force - check every pair
- FOR i in 1 .. (ST_NumPoints(ring)-2)
- LOOP
- FOR j in i .. (ST_NumPoints(ring)-1)
- LOOP
- d = ST_Distance(ST_PointN(ring,i),ST_PointN(ring,j));
- -- Check the distance and update if larger
- IF (d > dist) THEN
- dist = d;
- idx1 = i;
- idx2 = j;
- END IF;
- END LOOP;
- END LOOP;
-
- -- We now have the diameter of the convex hull. The following line returns it if desired.
- -- RETURN ST_MakeLine(ST_PointN(ring,idx1),ST_PointN(ring,idx2));
-
- -- Now for the Minimum Bounding Circle. Since we know the two points furthest from each
- -- other, the MBC must go through those two points. Start with those points as a diameter of a circle.
-
- -- The radius is half the distance between them and the center is midway between them
- radius = ST_Distance(ST_PointN(ring,idx1),ST_PointN(ring,idx2)) / 2.0;
- center = ST_LineInterpolatePoint(ST_MakeLine(ST_PointN(ring,idx1),ST_PointN(ring,idx2)),0.5);
-
- -- Loop through each vertex and check if the distance from the center to the point
- -- is greater than the current radius.
- FOR k in 1 .. (ST_NumPoints(ring)-1)
- LOOP
- IF(k <> idx1 and k <> idx2) THEN
- dist = ST_Distance(center,ST_PointN(ring,k));
- IF (dist > radius) THEN
- -- We have to expand the circle. The new circle must pass trhough
- -- three points - the two original diameters and this point.
-
- -- Draw a line from the first diameter to this point
- l1 = ST_Makeline(ST_PointN(ring,idx1),ST_PointN(ring,k));
- -- Compute the midpoint
- p1 = ST_LineInterpolatePoint(l1,0.5);
- -- Rotate the line 90 degrees around the midpoint (perpendicular bisector)
- l1 = ST_Rotate(l1,pi()/2,p1);
- -- Compute the azimuth of the bisector
- a1 = ST_Azimuth(ST_PointN(l1,1),ST_PointN(l1,2));
- -- Extend the line in each direction the new computed distance to insure they will intersect
- l1 = ST_AddPoint(l1,ST_Makepoint(ST_X(ST_PointN(l1,2))+sin(a1)*dist,ST_Y(ST_PointN(l1,2))+cos(a1)*dist),-1);
- l1 = ST_AddPoint(l1,ST_Makepoint(ST_X(ST_PointN(l1,1))-sin(a1)*dist,ST_Y(ST_PointN(l1,1))-cos(a1)*dist),0);
-
- -- Repeat for the line from the point to the other diameter point
- l2 = ST_Makeline(ST_PointN(ring,idx2),ST_PointN(ring,k));
- p2 = ST_LineInterpolatePoint(l2,0.5);
- l2 = ST_Rotate(l2,pi()/2,p2);
- a2 = ST_Azimuth(ST_PointN(l2,1),ST_PointN(l2,2));
- l2 = ST_AddPoint(l2,ST_Makepoint(ST_X(ST_PointN(l2,2))+sin(a2)*dist,ST_Y(ST_PointN(l2,2))+cos(a2)*dist),-1);
- l2 = ST_AddPoint(l2,ST_Makepoint(ST_X(ST_PointN(l2,1))-sin(a2)*dist,ST_Y(ST_PointN(l2,1))-cos(a2)*dist),0);
-
- -- The new center is the intersection of the two bisectors
- center = ST_Intersection(l1,l2);
- -- The new radius is the distance to any of the three points
- radius = ST_Distance(center,ST_PointN(ring,idx1));
- END IF;
- END IF;
- END LOOP;
- --DONE!! Return the MBC via the buffer command
- RETURN ST_Buffer(center,radius,segs_per_quarter);
-
- END;
-$BODY$
- LANGUAGE 'plpgsql' IMMUTABLE STRICT;
-
-
-- ST_ConcaveHull and Helper functions starts here --
-----------------------------------------------------------------------
-- Contributed by Regina Obe and Leo Hsu
Modified: trunk/regress/Makefile.in
===================================================================
--- trunk/regress/Makefile.in 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/regress/Makefile.in 2015-11-29 23:10:49 UTC (rev 14448)
@@ -95,6 +95,7 @@
long_xact \
lwgeom_regress \
measures \
+ minimum_bounding_circle \
operators \
out_geometry \
out_geography \
Modified: trunk/regress/tickets.sql
===================================================================
--- trunk/regress/tickets.sql 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/regress/tickets.sql 2015-11-29 23:10:49 UTC (rev 14448)
@@ -906,6 +906,12 @@
SELECT '#2956', st_astwkb(null,0) is null;
+-- #2996 --
+WITH
+ input AS (SELECT 'SRID=4326;POLYGON((26426 65078,26531 65242,26075 65136,26096 65427,26426 65078))'::geometry AS geom),
+ mbc AS (SELECT (mb).center, (mb).radius FROM (SELECT ST_MinimumBoundingRadius(geom) mb FROM input) sq)
+SELECT '#2996', radius = ST_Length(ST_LongestLine(geom, center)) FROM input, mbc;
+
SELECT '#3172', ST_AsText(ST_AddMeasure('LINESTRING(0 0,0 0)', 1, 2));
--SELECT '#3244a', ST_AsText(ST_3DClosestPoint('POINT(0 0 0)', 'POINT(0 0)'));
Modified: trunk/regress/tickets_expected
===================================================================
--- trunk/regress/tickets_expected 2015-11-28 08:31:14 UTC (rev 14447)
+++ trunk/regress/tickets_expected 2015-11-29 23:10:49 UTC (rev 14448)
@@ -275,6 +275,7 @@
#2788|f|Self-intersection|POINT(1 1)
#2870|Point[GS]
#2956|t
+#2996|t
#3172|LINESTRING M (0 0 1,0 0 2)
#3300|POLYGON((-71.7821 42.2622,-71.7821 42.9067,-71.029 42.9067,-71.029 42.2622,-71.7821 42.2622))
#3355|t
More information about the postgis-tickets
mailing list