[postgis-commits] svn - r3629 - trunk/lwgeom

[postgis-commits at postgis.refractions.net]

Author: robe
Date: 2009-02-03 06:21:25 -0800 (Tue, 03 Feb 2009)
New Revision: 3629

Modified:
   trunk/lwgeom/postgis.sql.in.c
Log:
First draft of ST_MinimumBoundingCircle contributed by Bruce Rindahl.  Changed to use named params and renamed function from mbc to ST_MinimumBoundingCircle.

Modified: trunk/lwgeom/postgis.sql.in.c
===================================================================
--- trunk/lwgeom/postgis.sql.in.c	2009-02-03 13:23:54 UTC (rev 3628)
+++ trunk/lwgeom/postgis.sql.in.c	2009-02-03 14:21:25 UTC (rev 3629)
@@ -6095,5 +6095,119 @@
 -- END
 ---------------------------------------------------------------
 
+
+---------------------------------------------------------------
+-- USER CONTRIUBUTED
+---------------------------------------------------------------
+
+-----------------------------------------------------------------------
+-- 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
+-----------------------------------------------------------------------
+CREATEFUNCTION ST_MinimumBoundingCircle(inputgeom geometry, segs_per_quarter integer)
+  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);
+	-- convert the hull perimeter to a linestring so we can manipulate individual points
+	ring = ST_ExteriorRing(hull);
+
+	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 MakeLine(PointN(ring,idx1),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_Line_interpolate_point(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_line_interpolate_point(l1,0.5);
+				-- Rotate the line 90 degrees around the midpoint (perpendicular bisector)
+				l1 = ST_Translate(ST_Rotate(ST_Translate(l1,-X(p1),-Y(p1)),pi()/2),X(p1),Y(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(X(ST_PointN(l1,2))+sin(a1)*dist,Y(ST_PointN(l1,2))+cos(a1)*dist),-1);
+				l1 = ST_AddPoint(l1,ST_Makepoint(X(ST_PointN(l1,1))-sin(a1)*dist,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_Line_interpolate_point(l2,0.5);
+				l2 = ST_Translate(ST_Rotate(ST_Translate(l2,-X(p2),-Y(p2)),pi()/2),X(p2),Y(p2));
+				a2 = ST_Azimuth(ST_PointN(l2,1),ST_PointN(l2,2));
+				l2 = ST_AddPoint(l2,ST_Makepoint(X(ST_PointN(l2,2))+sin(a2)*dist,Y(ST_PointN(l2,2))+cos(a2)*dist),-1);
+				l2 = ST_AddPoint(l2,ST_Makepoint(X(ST_PointN(l2,1))-sin(a2)*dist,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;
+  
+CREATE OR REPLACE FUNCTION ST_MinimumBoundingCircle(geometry)
+ RETURNS geometry AS
+'SELECT ST_MinimumBoundingCircle($1, 48)'
+ LANGUAGE 'sql' IMMUTABLE STRICT;
 COMMIT;