Initial work on ST_SplitGeometry. Split line by point implemented.

git-svn-id: http://svn.osgeo.org/postgis/trunk@5396 b70326c6-7e19-0410-871a-916f4a2858ee

postgis/Makefile.in

@@ -40,6 +40,7 @@ PG_OBJS=lwgeom_pg.o \
 	lwgeom_geos.o \
 	lwgeom_geos_prepared.o \
 	lwgeom_geos_clean.o \
+	lwgeom_geos_split.o \
 	lwgeom_export.o \
 	lwgeom_in_gml.o \
 	lwgeom_in_kml.o \

postgis/lwgeom_geos_split.c

@@ -0,0 +1,186 @@
+/**********************************************************************
+ * $Id: lwgeom_geos.c 5258 2010-02-17 21:02:49Z strk $
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ *
+ * Copyright 2009-2010 Sandro Santilli <strk@keybit.net>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************
+ *
+ * ST_SplitGeometry
+ *
+ * Split polygon by line, line by line, line by point.
+ * Returns components as a collection
+ *
+ * Author: Sandro Santilli <strk@keybit.net>
+ *
+ * Work done for Faunalia (http://www.faunalia.it) with fundings
+ * from Regione Toscana - Sistema Informativo per il Governo
+ * del Territorio e dell'Ambiente (RT-SIGTA).
+ *
+ * Thanks to the PostGIS community for sharing poly/line ideas [1]
+ *
+ * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
+ *
+ *
+ **********************************************************************/
+
+#include "lwgeom_geos.h"
+#include "funcapi.h"
+
+#include <string.h>
+#include <assert.h>
+
+/* #define POSTGIS_DEBUG_LEVEL 4 */
+
+/* Initializes GEOS internally */
+static LWGEOM* lwline_split_by_point(LWLINE* lwgeom_in, LWPOINT* blade_in);
+static LWGEOM*
+lwline_split_by_point(LWLINE* lwline_in, LWPOINT* blade_in)
+{
+	GEOSGeometry* line;
+	GEOSGeometry* point;
+	double loc, dist;
+	int intersects;
+	POINT2D pt;
+	POINTARRAY* pa1;
+	POINTARRAY* pa2;
+	LWGEOM** components;
+	LWCOLLECTION* out;
+
+	/* Possible outcomes:
+	 *
+	 *  1. The point is not on the line
+	 *      -> Return original geometry (cloned)
+	 *  2. The point is on the line
+	 *      -> Return a multiline with all elements
+	 */
+
+	getPoint2d_p(blade_in->point, 0, &pt);
+	loc = ptarray_locate_point(lwline_in->points, &pt, &dist);
+
+	/* lwnotice("Location: %g -- Distance: %g", loc, dist); */
+
+	if ( dist > 0 ) { /* TODO: accept a tolerance ? */
+		/* No intersection */
+		return (LWGEOM*)lwline_clone(lwline_in);
+	}
+
+	/* There is an intersection, let's get two substrings */
+
+	if ( loc == 0 || loc == 1 )
+	{
+		/* Intersection is on the boundary or outside */
+		return (LWGEOM*)lwline_clone(lwline_in);
+	}
+
+	pa1 = ptarray_substring(lwline_in->points, 0, loc);
+	pa2 = ptarray_substring(lwline_in->points, loc, 1);
+
+	/* TODO: check if either pa1 or pa2 are empty ? */
+
+	components = lwalloc(sizeof(LWGEOM*)*2);
+	components[0] = (LWGEOM*)lwline_construct(lwline_in->SRID, NULL, pa1);
+	components[1] = (LWGEOM*)lwline_construct(lwline_in->SRID, NULL, pa2);
+
+	out = lwcollection_construct(MULTILINETYPE, lwline_in->SRID,
+		NULL, 2, components);
+
+	/* That's all folks */
+
+	return (LWGEOM*)out;
+}
+
+static LWGEOM* lwline_split(LWLINE* lwgeom_in, LWGEOM* blade_in);
+static LWGEOM*
+lwline_split(LWLINE* lwline_in, LWGEOM* blade_in)
+{
+	switch (TYPE_GETTYPE(blade_in->type))
+	{
+	case POINTTYPE:
+		return lwline_split_by_point(lwline_in, (LWPOINT*)blade_in);
+
+	case LINETYPE:
+	default:
+		lwerror("Splitting a Line by a %s is unsupported",
+			lwgeom_typename(TYPE_GETTYPE(blade_in->type)));
+		return NULL;
+	}
+	return NULL;
+}
+
+static LWGEOM* lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in);
+static LWGEOM*
+lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in)
+{
+	switch (TYPE_GETTYPE(blade_in->type))
+	{
+	case LINETYPE:
+	default:
+		lwerror("Splitting a Polygon by a %s is unsupported",
+			lwgeom_typename(TYPE_GETTYPE(blade_in->type)));
+		return NULL;
+	}
+	return NULL;
+}
+
+static LWGEOM* lwgeom_split(LWGEOM* lwgeom_in, LWGEOM* blade_in);
+static LWGEOM*
+lwgeom_split(LWGEOM* lwgeom_in, LWGEOM* blade_in)
+{
+	switch (TYPE_GETTYPE(lwgeom_in->type))
+	{
+	case LINETYPE:
+		return lwline_split((LWLINE*)lwgeom_in, blade_in);
+
+	case POLYGONTYPE:
+		return lwpoly_split((LWPOLY*)lwgeom_in, blade_in);
+
+	default:
+		lwerror("Splitting of %d geometries is unsupported",
+			lwgeom_typename(TYPE_GETTYPE(lwgeom_in->type)));
+		return NULL;
+	}
+
+}
+
+
+Datum ST_SplitGeometry(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(ST_SplitGeometry);
+Datum ST_SplitGeometry(PG_FUNCTION_ARGS)
+{
+#if 0 && POSTGIS_GEOS_VERSION < 33
+	elog(ERROR, "You need GEOS-3.3.0 or up for ST_CleanGeometry");
+	PG_RETURN_NULL();
+#else /* POSTGIS_GEOS_VERSION >= 33 */
+
+	PG_LWGEOM *in, *blade_in, *out;
+	LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
+
+	in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	lwgeom_in = lwgeom_deserialize(SERIALIZED_FORM(in));
+
+	blade_in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
+	lwblade_in = lwgeom_deserialize(SERIALIZED_FORM(blade_in));
+
+	lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
+	if ( ! lwgeom_out ) {
+		PG_FREE_IF_COPY(in, 0);
+		PG_FREE_IF_COPY(blade_in, 1);
+		PG_RETURN_NULL();
+	}
+
+	out = pglwgeom_serialize(lwgeom_out);
+
+	PG_FREE_IF_COPY(in, 0);
+	PG_FREE_IF_COPY(blade_in, 1);
+
+	PG_RETURN_POINTER(out);
+
+#endif /* POSTGIS_GEOS_VERSION >= 33 */
+}
+

postgis/postgis.sql.in.c

@@ -4082,6 +4082,26 @@ CREATE OR REPLACE FUNCTION ST_CleanGeometry(geometry)
        LANGUAGE 'C' IMMUTABLE STRICT
        COST 100;
 
+--------------------------------------------------------------------------------
+-- ST_SplitGeometry 
+--------------------------------------------------------------------------------
+
+-- ST_SplitGeometry(in geometry, blade geometry)
+--
+-- Split a geometry in parts after cutting it with given blade.
+-- Returns a collection containing all parts.
+--
+-- Note that multi-part geometries will be returned exploded,
+-- no matter relation to blade.
+-- 
+-- Availability: 2.0.0
+--
+CREATE OR REPLACE FUNCTION ST_SplitGeometry(geometry, geometry)
+       RETURNS geometry
+       AS 'MODULE_PATHNAME', 'ST_SplitGeometry'
+       LANGUAGE 'C' IMMUTABLE STRICT
+       COST 100;
+
 --------------------------------------------------------------------------------
 -- Aggregates and their supporting functions
 --------------------------------------------------------------------------------

regress/Makefile.in

@@ -68,7 +68,8 @@ TESTS = \
 	dumppoints \
 	wmsservers \
 	tickets \
-	remove_repeated_points
+	remove_repeated_points \
+	split
 
 # Styled buffer only if GEOS >= 3.2
 ifeq ($(shell expr $(POSTGIS_GEOS_VERSION) ">=" 32),1)

regress/split.sql

@@ -0,0 +1,9 @@
+-- Point on line
+select '1',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(5 0)'));
+
+-- Point on line boundary
+select '2',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(10 0)'));
+
+-- Point off line
+select '3',st_astext(st_splitgeometry('LINESTRING(0 0, 10 0)', 'POINT(5 1)'));
+

regress/split_expected

@@ -0,0 +1,3 @@
+1|MULTILINESTRING((0 0,5 0),(5 0,10 0))
+2|LINESTRING(0 0,10 0)
+3|LINESTRING(0 0,10 0)