[postgis-users] Minimum Bounding Circle - was Roeck Test
Bruce Rindahl
rindahl at lrcwe.com
Thu Oct 30 08:27:09 PDT 2008
Robert
As mentioned by LEE and confirmed by me there are cases where this can
fail. The issue is in the first part of the code - finding the diameter
of the shape. This is not a trivial problem (Google "Diameter of a
Convex Polygon"). I can see why my code is failing when you look at the
edge cases but haven't yet found a solution. I am trying to avoid a
brute force approach (check the distances of every pair).
Bruce
Burgholzer,Robert wrote:
> 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
> rwburgholzer at deq.virginia.gov
> 804-698-4405
> Open Source Modeling Tools:
> http://sourceforge.net/projects/npsource/
>
> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net
> [mailto:postgis-users-bounces at postgis.refractions.net] 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
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
>
>
More information about the postgis-users
mailing list