Bruce, This is nice work! Certainly seems as if the problem is more complex than I had hoped!
r/b/ Robert W. Burgholzer Surface Water Modeler Office of Water Supply and Planning Virginia Department of Environmental Quality [EMAIL PROTECTED] 804-698-4405 Open Source Modeling Tools: http://sourceforge.net/projects/npsource/ -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Bruce Rindahl Sent: Wednesday, October 29, 2008 1:42 PM To: PostGIS Users Discussion Subject: [postgis-users] Minimum Bounding Circle - was Roeck Test The following is a procedure for a Minimum Bounding Circle of a geometry. It works with polygon geometry but it should also work with others. Anyone who is interested please test and comment. This was fun :-) Bruce Rindahl -- Function: mbc(geometry) -- DROP FUNCTION mbc(geometry); CREATE OR REPLACE FUNCTION mbc(geometry) 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 = convexHull($1); -- convert the hull perimeter to a linestring so we can manipulate individual points ring = exteriorRing(hull); -- Compute the diameter of the the convex hull -- The method is from Shamos (1978) -- The idea is to "roll" the polygon along the ground and find the maximum height of the bounding box -- Start with the side from the last point to the first point. p1 = PointN(ring,NumPoints(ring)-1); -- Set distance check to 0 dist = 0; FOR i in 1 .. (NumPoints(ring)-1) LOOP -- Compute the azimuth of the current side a1 = azimuth(p1,PointN(ring,i)); -- Rotate the polygon by 90 degrees + the azimuth this makes the side horizontal -- and compute the height of the bounding box d = ymax(box3d(rotate(ring,pi()/2+a1))) - ymin(box3d(rotate(ring,pi()/2+a1))); -- Check the distance and update if larger IF (d > dist) THEN dist = d; idx1 = i; END IF; -- Drag the current vertex along for the next side p1 = PointN(ring,i); END LOOP; -- The diameter (maximum distance between vertexes) -- is from either point idx1 or idx1-1. Compute the points at those indexes p1 = PointN(ring,idx1); IF (idx1 = 1) then p2 = PointN(ring,NumPoints(ring)-1); ELSE p2 = PointN(ring,idx1-1); END IF; -- Now find the vertex furthest from p1 or p2. The will be the -- other end of the diameter dist = 0; FOR j in 1 .. (NumPoints(ring)-1) LOOP IF (distance(PointN(ring,j),p1) > dist) THEN dist = distance(PointN(ring,j),p1); idx2 = j; END IF; IF (distance(PointN(ring,j),p2) > dist) THEN dist = distance(PointN(ring,j),p2); idx2 = j; END IF; END LOOP; -- idx2 now has the point index of the other end of the diameter -- Compute which is the starting point - p1 or p2 IF (distance(PointN(ring,idx2),p2) > distance(PointN(ring,idx2),p1)) THEN idx1 = idx1 - 1; END IF; IF (idx1 = 0) THEN idx1 = idx1 + NumPoints(ring) - 1; END IF; -- 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 = distance(PointN(ring,idx1),PointN(ring,idx2)) / 2.0; center = line_interpolate_point(MakeLine(PointN(ring,idx1),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 .. (NumPoints(ring)-1) LOOP IF(k <> idx1 and k <> idx2) THEN dist = distance(center,PointN(ring,k)); IF (dist > radius) THEN -- We have to expand the circle. The new circle must pass through -- three points - the two original diameters and this point. -- Draw a line from the first diameter to this point l1 = makeline(PointN(ring,idx1),PointN(ring,k)); -- Compute the midpoint p1 = line_interpolate_point(l1,0.5); -- Rotate the line 90 degrees around the midpoint (perpendicular bisector) l1 = translate(rotate(translate(l1,-X(p1),-Y(p1)),pi()/2),X(p1),Y(p1)); -- Compute the azimuth of the bisector a1 = azimuth(PointN(l1,1),PointN(l1,2)); -- Extend the line in each direction the new computed distance to insure they will intersect l1 = addPoint(l1,makepoint(X(PointN(l1,2))+sin(a1)*dist,Y(PointN(l1,2))+cos(a 1)*dist),-1); l1 = addPoint(l1,makepoint(X(PointN(l1,1))-sin(a1)*dist,Y(PointN(l1,1))-cos(a 1)*dist),0); -- Repeat for the line from the point to the other diameter point l2 = makeline(PointN(ring,idx2),PointN(ring,k)); p2 = line_interpolate_point(l2,0.5); l2 = translate(rotate(translate(l2,-X(p2),-Y(p2)),pi()/2),X(p2),Y(p2)); a2 = azimuth(PointN(l2,1),PointN(l2,2)); l2 = addPoint(l2,makepoint(X(PointN(l2,2))+sin(a2)*dist,Y(PointN(l2,2))+cos(a 2)*dist),-1); l2 = addPoint(l2,makepoint(X(PointN(l2,1))-sin(a2)*dist,Y(PointN(l2,1))-cos(a 2)*dist),0); -- The new center is the intersection of the two bisectors center = intersection(l1,l2); -- The new radius is the distance to any of the three points radius = distance(center,PointN(ring,idx1)); END IF; END IF; END LOOP; --DONE!! Return the MBC via the buffer command RETURN buffer(center,radius,48); END; $BODY$ LANGUAGE 'plpgsql' VOLATILE _______________________________________________ postgis-users mailing list [email protected] http://postgis.refractions.net/mailman/listinfo/postgis-users _______________________________________________ postgis-users mailing list [email protected] http://postgis.refractions.net/mailman/listinfo/postgis-users
