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
