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

[Lee Meilleur]

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

The district that is making it crash has an unusual shape:

http://www.gis.leg.mn/l2002/pdf/47A.pdf



-----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 12: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(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

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