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

[Lee Meilleur]

Bruce, it works great, thanks so much.  I'm sitting on a panel at an upcoming National Conference of State Legislatures meeting.  The panel is on Web Technologies for Redistricting and I'll be discussing PostgreSQL/PostGIS and MapServer as alternatives to some of the commercial software available.  Having a working MBC (Roeck test) is key.   Much appreciated.

Lee



-----Original Message-----
From: Bruce Rindahl [mailto:rindahl at lrcwe.com]
Sent: Thursday, October 30, 2008 11:24 AM
To: Lee Meilleur
Cc: 'PostGIS Users Discussion'
Subject: Re: [postgis-users] Minimum Bounding Circle - was Roeck Test



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