pggeodb.nancy.inra.fr/db_cefs - db_cefs on pggeodb.nancy.inra.fr

Function: st_minimumboundingcircle

Schema

public

Owner

postgres

Descriptions

args: geomA, num_segs_per_qt_circ=48 - Returns the smallest circle polygon that can fully contain a geometry. Default uses 48 segments per quarter circle.

Options

Definition

CREATE OR REPLACE FUNCTION public.st_minimumboundingcircle (
 inputgeom public.geometry,
 segs_per_quarter integer = 48
)
RETURNS public.geometry AS
$span$
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 = ST_ConvexHull(inputgeom);
--A point really has no MBC
IF ST_GeometryType(hull) = 'ST_Point' THEN
RETURN hull;
END IF;
-- convert the hull perimeter to a linestring so we can manipulate individual points
--If its already a linestring force it to a closed linestring
ring = CASE WHEN ST_GeometryType(hull) = 'ST_LineString' THEN ST_AddPoint(hull, ST_StartPoint(hull)) ELSE ST_ExteriorRing(hull) END;

dist = 0;
-- Brute Force - check every pair
FOR i in 1 .. (ST_NumPoints(ring)-2)
LOOP
FOR j in i .. (ST_NumPoints(ring)-1)
LOOP
d = ST_Distance(ST_PointN(ring,i),ST_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 ST_MakeLine(ST_PointN(ring,idx1),ST_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 = ST_Distance(ST_PointN(ring,idx1),ST_PointN(ring,idx2)) / 2.0;
center = ST_LineInterpolatePoint(ST_MakeLine(ST_PointN(ring,idx1),ST_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 .. (ST_NumPoints(ring)-1)
LOOP
IF(k <> idx1 and k <> idx2) THEN
dist = ST_Distance(center,ST_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 = ST_Makeline(ST_PointN(ring,idx1),ST_PointN(ring,k));
-- Compute the midpoint
p1 = ST_LineInterpolatePoint(l1,0.5);
-- Rotate the line 90 degrees around the midpoint (perpendicular bisector)
l1 = ST_Rotate(l1,pi()/2,p1);
--  Compute the azimuth of the bisector
a1 = ST_Azimuth(ST_PointN(l1,1),ST_PointN(l1,2));
--  Extend the line in each direction the new computed distance to insure they will intersect
l1 = ST_AddPoint(l1,ST_Makepoint(ST_X(ST_PointN(l1,2))+sin(a1)*dist,ST_Y(ST_PointN(l1,2))+cos(a1)*dist),-1);
l1 = ST_AddPoint(l1,ST_Makepoint(ST_X(ST_PointN(l1,1))-sin(a1)*dist,ST_Y(ST_PointN(l1,1))-cos(a1)*dist),0);

-- Repeat for the line from the point to the other diameter point
l2 = ST_Makeline(ST_PointN(ring,idx2),ST_PointN(ring,k));
p2 = ST_LineInterpolatePoint(l2,0.5);
l2 = ST_Rotate(l2,pi()/2,p2);
a2 = ST_Azimuth(ST_PointN(l2,1),ST_PointN(l2,2));
l2 = ST_AddPoint(l2,ST_Makepoint(ST_X(ST_PointN(l2,2))+sin(a2)*dist,ST_Y(ST_PointN(l2,2))+cos(a2)*dist),-1);
l2 = ST_AddPoint(l2,ST_Makepoint(ST_X(ST_PointN(l2,1))-sin(a2)*dist,ST_Y(ST_PointN(l2,1))-cos(a2)*dist),0);

-- The new center is the intersection of the two bisectors
center = ST_Intersection(l1,l2);
-- The new radius is the distance to any of the three points
radius = ST_Distance(center,ST_PointN(ring,idx1));
END IF;
END IF;
END LOOP;
--DONE!!  Return the MBC via the buffer command
RETURN ST_Buffer(center,radius,segs_per_quarter);

END;
$span$
LANGUAGE 'plpgsql'
IMMUTABLE
RETURNS NULL ON NULL INPUT
SECURITY INVOKER
COST 100;

COMMENT ON FUNCTION public.st_minimumboundingcircle(inputgeom public.geometry, segs_per_quarter integer)
IS 'args: geomA, num_segs_per_qt_circ=48 - Returns the smallest circle polygon that can fully contain a geometry. Default uses 48 segments per quarter circle.';

This file was generated with SQL Manager for PostgreSQL (www.pgsqlmanager.com) at 13/03/2014 13:23