# [postgis-users] Minimum Bounding Circle - was Roeck Test

Burgholzer,Robert rwburgholzer at deq.virginia.gov
Thu Oct 30 07:28:32 PDT 2008

```Bruce,
This is nice work!  Certainly seems as if the problem is more complex

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;
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
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));
-- 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 =
1)*dist),-1);
l1 =
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 =
2)*dist),-1);
l2 =
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
END IF;
END IF;
END LOOP;
--DONE!!  Return the MBC via the buffer command

END;
\$BODY\$
LANGUAGE 'plpgsql' VOLATILE

_______________________________________________
postgis-users mailing list
postgis-users at postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

```