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

[Bruce Rindahl]

Lee Meilleur wrote:
> Thanks Bruce, it works well.  When I compare it with our commercial software it's within a couple of hundredths of the Roeck test, but not in any predictable fashion.  Some are .02 over the Roeck, some are .01 under,  the majority are right on, but there were several that were  1 to 3 tenths off, fairly significant.  I have no way of knowing how the commercial software computes the Roeck test.  There was one district that created an error, I had to run the query on each of the districts to figure it out.  Here's the message:
>
> Error: getPoint2d_p:  point offset out of range
>
>   
Lee and anyone else interested.

There are two issues to be resolved here.  The first is the one that I 
mentioned previously - the computation of the diameter of a convex 
polygon efficiently is a nasty problem.  I am looking into correcting 
the original code.  Below is a revised version that will work by brute 
force - it checks every pair of points.  If the convex polygon doesn't 
have too many points it works well.  Increasing the number of points 
increases the computation time exponentially.  You can test this by the 
following simple query:

SELECT mbc(buffer(MakePoint(0,0),100,8)) from mytable

Buffering a point actually creates a polygon - not a circle.  The number 
of points in the polygon is 4 times the last parameter (8) in the above 
query. Eight is the default.  Watch the computation time increase as you 
increase the value of 8 in the query.  However Lee's data set and mine 
run very fast using the Brute Force approach.

The second issue is the different results of the area of the circle.  I 
think this is again the result of a circle is actually a polygon.  In 
the code the function returns a buffer around the computed center with 
the computed radius.  I set the number of points to 48 * 4 - see the 
RETURN line.  This is probably close for small areas but it may need to 
be increased if you need a polygon result.  For an exact calculation, 
the radius could be returned (or an area) if desired.

Bruce

-- 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);

    dist = 0;
    -- Brute Force - check every pair
    FOR i in 1 .. (NumPoints(ring)-2)
        LOOP
            FOR j in i .. (NumPoints(ring)-1)
                LOOP
                d = distance(PointN(ring,i),PointN(ring,j));
                -- Check the distance and update if larger
                IF (d > dist) THEN
                    dist = d;
                    idx1 = i;
                    idx2 = j;
                END IF;
            END LOOP;
        END LOOP;

    -- 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 trhough
                -- 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