[prev in list] [next in list] [prev in thread] [next in thread] 

List:       postgis-users
Subject:    Re: [postgis-users] Cutting a Polygon into Pieces with maximum
From:       Markus Schaber <schabi () logix-tt ! com>
Date:       2008-06-24 9:18:45
Message-ID: 20080624111845.4fe221f5 () kingfisher ! sec ! intern ! logix-tt ! com
[Download RAW message or body]

Hi, all,

I now have an basic implementation of this functionality in plpgsql,
which I post here for all those googlers who might have the same
problem.

I'm not shure whether it is in shape for committing, yet.


Regards,
Markus

-- 
Markus Schaber | Logical Tracking&Tracing International AG
Dipl. Inf.     | Software Development GIS

Fight against software patents in Europe! www.ffii.org
www.nosoftwarepatents.org

["prep_2_addPolySplitter.sql" (text/x-sql)]

--
-- Function and helpers to split a Polygon into multiple smaller polygons (fragments).
-- © 2008 Markus Schaber <schabi@logix-tt.com>
-- Licensed under GPL V2 or later.

-- Helper function to get a point from a linestring with wraparound at the end.
CREATE OR REPLACE FUNCTION pointNWrap(in_line geometry, -- LINESTRING 
					index int4) -- Index based on 1
RETURNS geometry -- POINT
LANGUAGE SQL
IMMUTABLE
AS $$
	SELECT PointN($1, ($2-1)%nPoints($1) + 1);
$$;

-- Helper function to get substring from a linestring with wraparound at the end.
-- This could be speed up with a native implementation like the pgis_GeometryRangeN()
-- discussed on the PostGIS list.
CREATE OR REPLACE FUNCTION splitLineStringWrap(in_line geometry, -- LINESTRING
						start int4, -- Indices based on 1 
						len int4) -- Number of points, will be wrapped / expanded.
RETURNS geometry -- LINESTRING
LANGUAGE SQL
IMMUTABLE
AS $$
SELECT makeLine(point) FROM (
	SELECT pointNWrap($1, generate_series($2, $2+$3-1)) AS point) 
AS foo;
$$;

-- This function splits a polygon into partial polygons which have at
-- maximum maxPoints points in their outer ring.
-- This algorithm does not guarantee the optimal result, it may deliver
-- more (and smaller) fragments than necessary.
-- Currently, all inner Rings are ignored (dropped), this could clearly
-- be enhanced in a future version. It is even possible to make an
-- implementation that guarantees to eliminate all inner Rings by
-- cutting appropriately.
-- The current implementation also contains some debug output which
-- was commented out for speed reasons.
CREATE OR REPLACE FUNCTION _splitPolyMaxVertex(in_poly geometry, maxPoints int4)
RETURNS SETOF geometry
LANGUAGE plpgsql
IMMUTABLE
AS $$
DECLARE
	currGap int4;
	currIndex int4;
	closedRing geometry;
	openRing geometry;
	restPoly geometry;
	
	cutStart geometry;
	cutEnd geometry;
	cutLine geometry;
	
	splitPart geometry;
	splitPoly geometry;
BEGIN
--	RAISE NOTICE ' +++ next polygon +++ ';
	IF geometryType(in_poly) != 'POLYGON' THEN
		RAISE WARNING 'Dropping bad geometry type %', geometryType(in_poly);
		return;
	ELSIF not isvalid(in_poly) THEN
		RAISE WARNING 'Dropping invalid geometry';
		return;
	ELSIF maxPoints <3 THEN
		RAISE EXCEPTION 'maxPoints=%, but must be 4 or greater!', maxPoints;
	ELSIF nPoints(exteriorRing(in_poly))<=maxPoints THEN
		-- polygon is small enough
		return next in_poly;
		return;
	ELSIF NumInteriorRings(in_poly)>0 THEN
		RAISE NOTICE 'Ignoring % inner rings', NumInteriorRings(in_poly);
	END IF;
	
	closedRing := exteriorring(in_poly);
	restPoly := makePolygon(closedRing);
	openRing := removePoint(closedRing, nPoints(closedRing)-1);

	currIndex := 1;
	currGap := maxPoints-1;
	
--	RAISE DEBUG 'Exterior ring before outer loop: % %', nPoints(closedRing), asText(closedRing);
	WHILE currGap > 0 LOOP
--		RAISE NOTICE ' --- next try --- ';
		cutStart := pointN(openRing, currIndex);
		cutEnd := pointNWrap(openRing, currIndex + currGap-1);		
		cutLine := makeLine(cutStart, cutEnd);
		
		IF contains(restPoly, cutLine) THEN
--			RAISE NOTICE ' --- next cut --- ';
			splitPart := splitLineStringWrap(openRing, currIndex, currGap);
			
--			IF startPoint(splitPart)!=cutStart THEN
--				RAISE WARNING 'cut start point broken. % %', asText(cutLine), asText(splitPart);
--			END IF;

--			IF endPoint(splitPart) != cutend THEN
--				RAISE WARNING 'cut end point broken. % %', asText(cutLine), asText(splitPart);
--			END IF;
			
			splitPoly := makePolygon(addPoint(splitPart, startPoint(splitPart)));
			
--			IF not isvalid(splitPoly) THEN
--				RAISE WARNING 'splitPoly invalid.';
--			END IF;
			return next splitPoly;
			
			openRing := splitLineStringWrap(openRing, currIndex+currGap-1, nPoints(openRing)-currGap+2);

--			IF startPoint(openRing)!=cutEnd THEN
--				RAISE WARNING 'rest start point broken. % %', asText(cutLine), asText(openRing);
--			END IF;
			
--			IF endPoint(openRing) != cutStart THEN
--				RAISE WARNING 'rest end point broken. % %', asText(cutLine), asText(openRing);
--			END IF;
			
			closedRing := addPoint(openRing, startPoint(openRing));
			restPoly := makePolygon(closedRing);
			
--			IF not isvalid(restpoly) THEN
--				RAISE WARNING 'restPoly invalid.';
--			END IF;

			IF nPoints(closedRing) <= maxPoints THEN
				return next restPoly;
				return;
			END IF;

			currIndex := 1;
			currGap := maxPoints -1; -- little slowdown, but potentially delivers better results.
		ELSIF currIndex = nPoints(openRing) THEN
			currIndex := 1;
			currGap := currGap-1;
		ELSE
			currIndex := currIndex + 1;					 
		END IF;
	END LOOP;
	RAISE EXCEPTION 'ASSERTION failed - could not cut polygon into triangles!';	      
END;
$$;

-- SQL language wrapper for _splitPolyMaxVertex.
-- Due to some historical artifacts in the plpgsql implementation, there are some
-- limitations in what contexts set returning functions can be called. Those do not
-- apply to SQL functions, so creating an SQL wrapper around the plpgsql function
-- solves the problem.
CREATE OR REPLACE FUNCTION splitPolyMaxVertex(in_poly geometry, maxPoints int4)
RETURNS SETOF geometry
LANGUAGE sql 
IMMUTABLE
AS $$
	SELECT * FROM _splitPolyMaxVertex($1, $2);
$$;


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


[prev in list] [next in list] [prev in thread] [next in thread]