-- Contributed by Regina Obe and Leo Hsu\r
-- Availability: 2.0.0\r
-----------------------------------------------------------------------\r
-CREATE OR REPLACE FUNCTION _ST_ConcaveHull(param_inputgeom geometry)\r
+CREATE OR REPLACE FUNCTION _st_concavehull(param_inputgeom geometry)\r
RETURNS geometry AS\r
$$\r
DECLARE \r
vexhull GEOMETRY;\r
var_resultgeom geometry;\r
+ var_inputgeom geometry;\r
vexring GEOMETRY;\r
cavering GEOMETRY;\r
cavept geometry[];\r
seglength double precision;\r
var_tempgeom geometry;\r
+ scale_factor integer := 1;\r
i integer;\r
\r
BEGIN\r
\r
-- First compute the ConvexHull of the geometry\r
vexhull := ST_ConvexHull(param_inputgeom);\r
+ var_inputgeom := param_inputgeom;\r
--A point really has no concave hull\r
IF ST_GeometryType(vexhull) = 'ST_Point' OR ST_GeometryType(vexHull) = 'ST_LineString' THEN\r
RETURN vexhull;\r
\r
-- convert the hull perimeter to a linestring so we can manipulate individual points\r
vexring := CASE WHEN ST_GeometryType(vexhull) = 'ST_LineString' THEN vexhull ELSE ST_ExteriorRing(vexhull) END;\r
+ IF abs(ST_X(ST_PointN(vexring,1))) < 1 THEN --scale the geometry to prevent stupid precision errors - not sure it works so make low for now\r
+ scale_factor := 100;\r
+ vexring := ST_Scale(vexring, scale_factor,scale_factor);\r
+ var_inputgeom := ST_Scale(var_inputgeom, scale_factor, scale_factor);\r
+ --RAISE NOTICE 'Scaling';\r
+ END IF;\r
seglength := ST_Length(vexring)/least(ST_NPoints(vexring)*2,1000) ;\r
\r
vexring := ST_Segmentize(vexring, seglength);\r
ARRAY(\r
\r
SELECT \r
- ST_ClosestPoint(param_inputgeom, pt ) As the_geom\r
+ ST_ClosestPoint(var_inputgeom, pt ) As the_geom\r
FROM (\r
SELECT ST_PointN(vexring, n ) As pt, n\r
FROM \r
)\r
)\r
; \r
+ \r
\r
- var_resultgeom := ST_MakePolygon(ST_MakeLine(geom)) \r
+ var_resultgeom := ST_MakeLine(geom) \r
FROM ST_Dump(cavering) As foo;\r
+\r
+ IF ST_IsSimple(var_resultgeom) THEN\r
+ var_resultgeom := ST_MakePolygon(var_resultgeom);\r
+ --RAISE NOTICE 'is Simple: %', var_resultgeom;\r
+ ELSE /** will not result in a valid polygon -- just return convex hull **/\r
+ --RAISE NOTICE 'is not Simple: %', var_resultgeom;\r
+ var_resultgeom := ST_ConvexHull(var_resultgeom);\r
+ END IF;\r
\r
- IF NOT ST_IsValid(var_resultgeom) THEN\r
- --RAISE NOTICE '_ST_Concavehull invalid %', ST_AsText(var_resultgeom);\r
- var_tempgeom := ST_BuildArea(var_resultgeom); -- try to make valid\r
- IF NOT ST_IsValid(var_tempgeom) THEN\r
- var_resultgeom := ST_Buffer(var_resultgeom,ST_Length(cavering)/1000, 'quad_segs=3'); -- try to make valid\r
- END IF;\r
- --if still invalid or doens't contain the geometry just return convex hull\r
- IF NOT ST_IsValid(var_resultgeom) or ST_GeometryType(var_resultgeom) <> 'ST_Polygon' THEN\r
- var_resultgeom := ST_ConvexHull(param_inputgeom);\r
- ELSIF ST_GeometryType(param_inputgeom) ILIKE '%Geometry%' THEN\r
- IF EXISTS(SELECT geom FROM ST_Dump(param_inputgeom) WHERE NOT ST_Covers(var_resultgeom,geom) ) THEN \r
- --we have to explode inputgeom since geos doesn't support geometrycollections for containment check\r
- var_resultgeom := ST_ConvexHull(param_inputgeom);\r
- END IF;\r
- ELSIF NOT ST_Contains(var_resultgeom, param_inputgeom) THEN\r
- var_resultgeom := ST_ConvexHull(param_inputgeom);\r
- END IF;\r
+ IF scale_factor > 1 THEN -- scale the result back\r
+ var_resultgeom := ST_Scale(var_resultgeom, 1/scale_factor, 1/scale_factor);\r
END IF;\r
RETURN var_resultgeom;\r
\r
END;\r
$$\r
- LANGUAGE 'plpgsql' IMMUTABLE STRICT\r
- COST 100;\r
+ LANGUAGE plpgsql IMMUTABLE STRICT;\r
\r
CREATE OR REPLACE FUNCTION ST_ConcaveHull(param_geom geometry, param_pctconvex float, param_allow_holes boolean) RETURNS geometry AS\r
$$\r
--- /dev/null
+-- $Id$\r
+-- Tests to confirm the concave hull area is <= convex hull and \r
+-- covers the original geometry (can't use covers because always gives topo errors with 3.3\r
+SELECT \r
+ 'ST_ConcaveHull MultiPolygon 0.95', ST_Area(ST_Intersection(geom,ST_ConcaveHull(\r
+ geom, 0.95) )) = ST_Area(geom) As encloses_geom, \r
+ (ST_Area(ST_ConvexHull(geom)) \r
+ - ST_Area(ST_ConcaveHull(geom, 0.95))) < (0.95 * ST_Area(ST_ConvexHull(geom) ) ) As reached_target\r
+FROM ST_Union(ST_GeomFromText('POLYGON((175 150, 20 40, \r
+ 50 60, 125 100, 175 150))'),\r
+ ST_Buffer(ST_GeomFromText('POINT(110 170)'), 20)\r
+ ) As geom;\r
+ \r
+SELECT \r
+ 'ST_ConcaveHull Lines 0.80', ST_Intersection(geom,ST_ConcaveHull(\r
+ geom, 0.80) ) = geom As encloses_geom, \r
+ (ST_Area(ST_ConvexHull(geom)) \r
+ - ST_Area(ST_ConcaveHull(geom, 0.80))) < (0.80 * ST_Area(ST_ConvexHull(geom) ) ) As reached_target\r
+\r
+FROM ST_GeomFromText('MULTILINESTRING((106 164,30 112,74 70,82 112,130 94,\r
+ 130 62,122 40,156 32,162 76,172 88),\r
+(132 178,134 148,128 136,96 128,132 108,150 130,\r
+170 142,174 110,156 96,158 90,158 88),\r
+(22 64,66 28,94 38,94 68,114 76,112 30,\r
+132 10,168 18,178 34,186 52,184 74,190 100,\r
+190 122,182 148,178 170,176 184,156 164,146 178,\r
+132 186,92 182,56 158,36 150,62 150,76 128,88 118))') As geom;\r
+\r
+-- test holes vs. no holes - holes should still enclose but have smaller area than no holes --\r
+SELECT \r
+ 'ST_ConcaveHull Lines 0.80 holes', ST_Intersection(geom,ST_ConcaveHull(\r
+ geom, 0.80, true) ) = geom As encloses_geom, \r
+ ST_Area(ST_ConcaveHull(geom, 0.80, true)) < ST_Area(ST_ConcaveHull(geom, 0.80)) As reached_target\r
+\r
+FROM ST_GeomFromText('MULTILINESTRING((106 164,30 112,74 70,82 112,130 94,\r
+ 130 62,122 40,156 32,162 76,172 88),\r
+(132 178,134 148,128 136,96 128,132 108,150 130,\r
+170 142,174 110,156 96,158 90,158 88),\r
+(22 64,66 28,94 38,94 68,114 76,112 30,\r
+132 10,168 18,178 34,186 52,184 74,190 100,\r
+190 122,182 148,178 170,176 184,156 164,146 178,\r
+132 186,92 182,56 158,36 150,62 150,76 128,88 118))') As geom;\r
projects / postgis / commitdiff