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

[Bruce Rindahl]

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(a1)*dist),-1);
                l1 = 
addPoint(l1,makepoint(X(PointN(l1,1))-sin(a1)*dist,Y(PointN(l1,1))-cos(a1)*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(a2)*dist),-1);
                l2 = 
addPoint(l2,makepoint(X(PointN(l2,1))-sin(a2)*dist,Y(PointN(l2,1))-cos(a2)*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