-- Function: mbc(geometry, integer) -- DROP FUNCTION mbc(geometry, integer); CREATE OR REPLACE FUNCTION mbc(geometry, 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(\$1); -- 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,\$2); END; \$BODY\$ LANGUAGE 'plpgsql' VOLATILE COST 100; ALTER FUNCTION mbc(geometry, integer) OWNER TO postgres; -- Function: mbc(geometry) -- DROP FUNCTION mbc(geometry); CREATE OR REPLACE FUNCTION mbc(geometry) RETURNS geometry AS 'SELECT mbc(\$1, 48)' LANGUAGE 'sql' IMMUTABLE STRICT COST 100; ALTER FUNCTION mbc(geometry) OWNER TO postgres;