[postgis-commits] svn - r3547 - in trunk: . lwgeom
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Wed Jan 21 15:19:50 PST 2009
Author: pramsey
Date: 2009-01-21 15:19:50 -0800 (Wed, 21 Jan 2009)
New Revision: 3547
Modified:
trunk/configure.ac
trunk/lwgeom/lwgeom_geos.c
trunk/lwgeom/postgis.sql.in.c
Log:
Add support for fast unions, with cascaded union. Currently for testing, in the ST_Union_Fast() agggregate. Requires GEOS SVN r2252 or higher.
Modified: trunk/configure.ac
===================================================================
--- trunk/configure.ac 2009-01-21 21:55:13 UTC (rev 3546)
+++ trunk/configure.ac 2009-01-21 23:19:50 UTC (rev 3547)
@@ -447,7 +447,7 @@
LIBS="$LIBS_SAVE"
dnl ===========================================================================
-dnl Detect GLib GTK for GUI
+dnl Detect GTK+2.0 for GUI
dnl ===========================================================================
AC_ARG_WITH([gui],
Modified: trunk/lwgeom/lwgeom_geos.c
===================================================================
--- trunk/lwgeom/lwgeom_geos.c 2009-01-21 21:55:13 UTC (rev 3546)
+++ trunk/lwgeom/lwgeom_geos.c 2009-01-21 23:19:50 UTC (rev 3547)
@@ -47,6 +47,7 @@
Datum symdifference(PG_FUNCTION_ARGS);
Datum geomunion(PG_FUNCTION_ARGS);
Datum unite_garray(PG_FUNCTION_ARGS);
+Datum unite_garray_fast(PG_FUNCTION_ARGS);
Datum issimple(PG_FUNCTION_ARGS);
Datum isring(PG_FUNCTION_ARGS);
Datum geomequals(PG_FUNCTION_ARGS);
@@ -81,7 +82,6 @@
PG_RETURN_POINTER(result);
}
-#ifndef UNITE_USING_BUFFER
/*
* This is the final function for GeomUnion
* aggregate. Will have as input an array of Geometries.
@@ -89,8 +89,8 @@
* versions of them and return PGIS-converted version back.
* Changing combination order *might* speed up performance.
*/
-PG_FUNCTION_INFO_V1(unite_garray);
-Datum unite_garray(PG_FUNCTION_ARGS)
+PG_FUNCTION_INFO_V1(unite_garray_fast);
+Datum unite_garray_fast(PG_FUNCTION_ARGS)
{
Datum datum;
ArrayType *array;
@@ -99,10 +99,13 @@
PG_LWGEOM *result, *pgis_geom;
GEOSGeom g1, g2, geos_result=NULL;
int SRID=-1;
- size_t offset;
+ size_t offset = 0;
#if POSTGIS_DEBUG_LEVEL > 0
static int call=1;
#endif
+#if POSTGIS_GEOS_VERSION >= 31
+ int allpolys=1;
+#endif
#if POSTGIS_DEBUG_LEVEL > 0
call++;
@@ -128,96 +131,179 @@
/* Ok, we really need geos now ;) */
initGEOS(lwnotice, lwnotice);
+#if POSTGIS_GEOS_VERSION >= 31
+
+ /*
+ ** First, see if all our elements are POLYGON/MULTIPOLYGON
+ ** If they are, we can use UnionCascaded for faster results.
+ */
offset = 0;
- for (i=0; i<nelems; i++)
+ for ( i = 0; i < nelems; i++ )
{
- PG_LWGEOM *geom = (PG_LWGEOM *)(ARR_DATA_PTR(array)+offset);
- offset += INTALIGN(VARSIZE(geom));
+ PG_LWGEOM *pggeom = (PG_LWGEOM *)(ARR_DATA_PTR(array)+offset);
+ int pgtype = TYPE_GETTYPE(pggeom->type);
+ offset += INTALIGN(VARSIZE(pggeom));
+ if( ! i ) /* Initialize SRID */
+ {
+ SRID = pglwgeom_getSRID(pggeom);
+ if ( TYPE_HASZ(pggeom->type) ) is3d = 1;
+ }
+ if ( pgtype != POLYGONTYPE && pgtype != MULTIPOLYGONTYPE )
+ {
+ allpolys = 0;
+ break;
+ }
+ }
- pgis_geom = geom;
+ if ( allpolys )
+ {
+ /*
+ ** Everything is polygonal, so let's run the cascaded polygon union!
+ */
+ int geoms_size = nelems;
+ int curgeom = 0;
+ GEOSGeom *geoms = NULL;
+ geoms = palloc( sizeof(GEOSGeom) * geoms_size );
+ /*
+ ** We need to convert the array of PG_LWGEOM into a GEOS MultiPolygon.
+ ** First make an array of GEOS Polygons.
+ */
+ offset = 0;
+ for( i = 0; i < nelems; i++ )
+ {
+ PG_LWGEOM *pggeom = (PG_LWGEOM *)(ARR_DATA_PTR(array)+offset);
+ int pgtype = TYPE_GETTYPE(pggeom->type);
+ offset += INTALIGN(VARSIZE(pggeom));
+ if( pgtype == POLYGONTYPE )
+ {
+ if( curgeom == geoms_size )
+ {
+ geoms_size *= 2;
+ geoms = repalloc( geoms, sizeof(GEOSGeom) * geoms_size );
+ }
+ geoms[curgeom] = POSTGIS2GEOS(pggeom);
+ curgeom++;
+ }
+ if( pgtype == MULTIPOLYGONTYPE )
+ {
+ int j = 0;
+ LWGEOM_INSPECTED *lwgeom = lwgeom_inspect(SERIALIZED_FORM(pggeom));
+ for ( j = 0; j < lwgeom->ngeometries; j++ )
+ {
+ LWPOLY *lwpoly = NULL;
+ int k = 0;
+ if( curgeom == geoms_size )
+ {
+ geoms_size *= 2;
+ geoms = repalloc( geoms, sizeof(GEOSGeom) * geoms_size );
+ }
+ /* This builds a LWPOLY on top of the serialized form */
+ lwpoly = lwgeom_getpoly_inspected(lwgeom, j);
+ geoms[curgeom] = LWGEOM2GEOS(lwpoly_as_lwgeom(lwpoly));
+ /* We delicately free the LWPOLY and POINTARRAY structs, leaving the serialized form below untouched. */
+ for( k = 0; k < lwpoly->nrings; k++ )
+ {
+ lwfree(lwpoly->rings[k]);
+ }
+ lwpoly_release(lwpoly);
+ curgeom++;
+ }
+ }
- POSTGIS_DEBUGF(3, "geom %d @ %p", i, geom);
+ }
+ /*
+ ** Take our GEOS Polygons and turn them into a GEOS MultiPolygon,
+ ** then pass that into cascaded union.
+ */
+ g1 = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, curgeom);
+ /* TODO protect against null return */
+ g2 = GEOSUnionCascaded(g1);
+ /* TODO protect against null return */
+ GEOSSetSRID(g2, SRID);
+ result = GEOS2POSTGIS(g2, is3d);
+ /* Clean up the mess. */
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(g2);
+ }
+ else
+ {
+#endif /* POSTGIS_GEOS_VERSION >= 31 */
+ /*
+ ** Heterogeneous result, let's slog through this one union at a time.
+ */
+ offset = 0;
+ for (i=0; i<nelems; i++)
+ {
+ PG_LWGEOM *geom = (PG_LWGEOM *)(ARR_DATA_PTR(array)+offset);
+ offset += INTALIGN(VARSIZE(geom));
- /* Check is3d flag */
- if ( TYPE_HASZ(geom->type) ) is3d = 1;
+ pgis_geom = geom;
- /* Check SRID homogeneity and initialize geos result */
- if ( ! i )
- {
- geos_result = POSTGIS2GEOS(geom);
- SRID = pglwgeom_getSRID(geom);
+ POSTGIS_DEBUGF(3, "geom %d @ %p", i, geom);
- POSTGIS_DEBUGF(3, "first geom is a %s", lwgeom_typename(TYPE_GETTYPE(geom->type)));
+ /* Check is3d flag */
+ if ( TYPE_HASZ(geom->type) ) is3d = 1;
- continue;
- }
- else
- {
- errorIfSRIDMismatch(SRID, pglwgeom_getSRID(geom));
- }
+ /* Check SRID homogeneity and initialize geos result */
+ if ( ! i )
+ {
+ geos_result = POSTGIS2GEOS(geom);
+ SRID = pglwgeom_getSRID(geom);
+ POSTGIS_DEBUGF(3, "first geom is a %s", lwgeom_typename(TYPE_GETTYPE(geom->type)));
+ continue;
+ }
+ else
+ {
+ errorIfSRIDMismatch(SRID, pglwgeom_getSRID(geom));
+ }
- g1 = POSTGIS2GEOS(pgis_geom);
+ g1 = POSTGIS2GEOS(pgis_geom);
- POSTGIS_DEBUGF(3, "unite_garray(%d): adding geom %d to union (%s)",
- call, i, lwgeom_typename(TYPE_GETTYPE(geom->type)));
+ POSTGIS_DEBUGF(3, "unite_garray(%d): adding geom %d to union (%s)",
+ call, i, lwgeom_typename(TYPE_GETTYPE(geom->type)));
- g2 = GEOSUnion(g1,geos_result);
- if ( g2 == NULL )
- {
+ g2 = GEOSUnion(g1,geos_result);
+ if ( g2 == NULL )
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(geos_result);
+ elog(ERROR,"GEOS union() threw an error!");
+ }
GEOSGeom_destroy(g1);
GEOSGeom_destroy(geos_result);
- elog(ERROR,"GEOS union() threw an error!");
+ geos_result = g2;
}
- GEOSGeom_destroy(g1);
+
+ GEOSSetSRID(geos_result, SRID);
+ result = GEOS2POSTGIS(geos_result, is3d);
GEOSGeom_destroy(geos_result);
- geos_result = g2;
+
+#if POSTGIS_GEOS_VERSION >= 31
}
-
- GEOSSetSRID(geos_result, SRID);
- result = GEOS2POSTGIS(geos_result, is3d);
- GEOSGeom_destroy(geos_result);
+#endif
+
if ( result == NULL )
{
elog(ERROR, "GEOS2POSTGIS returned an error");
PG_RETURN_NULL(); /* never get here */
}
- /* compressType(result); */
-
PG_RETURN_POINTER(result);
}
-#else /* def UNITE_USING_BUFFER */
-
-/*
-* This is the final function for GeomUnion
-* aggregate. Will have as input an array of Geometries.
-* Builds a GEOMETRYCOLLECTION from input and call
-* GEOSBuffer(collection, 0) on the GEOS-converted
-* versions of it. Returns PGIS-converted version back.
-*/
PG_FUNCTION_INFO_V1(unite_garray);
Datum unite_garray(PG_FUNCTION_ARGS)
{
Datum datum;
ArrayType *array;
int is3d = 0;
- int nelems, i, ngeoms, npoints;
- PG_LWGEOM *result=NULL;
- GEOSGeom *geoms, collection;
- GEOSGeom g1, geos_result=NULL;
+ int nelems, i;
+ PG_LWGEOM *result, *pgis_geom;
+ GEOSGeom g1, g2, geos_result=NULL;
int SRID=-1;
- size_t offset;
-#if POSTGIS_DEBUG_LEVEL > 0
- static int call=1;
-#endif
+ size_t offset = 0;
-#if POSTGIS_DEBUG_LEVEL >= 2
- call++;
- POSTGIS_DEBUGF(2, "GEOS buffer union (call %d)", call);
-#endif
-
datum = PG_GETARG_DATUM(0);
/* Null array, null geometry (should be empty?) */
@@ -227,106 +313,59 @@
nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
- POSTGIS_DEBUGF(3, "unite_garray: number of elements: %d", nelems);
-
if ( nelems == 0 ) PG_RETURN_NULL();
/* One-element union is the element itself */
if ( nelems == 1 ) PG_RETURN_POINTER((PG_LWGEOM *)(ARR_DATA_PTR(array)));
- geoms = lwalloc(sizeof(GEOSGeom)*nelems);
-
- /* We need geos here */
+ /* Ok, we really need geos now ;) */
initGEOS(lwnotice, lwnotice);
offset = 0;
- i=0;
- ngeoms = 0;
- npoints=0;
-
- POSTGIS_DEBUGF(3, "Nelems %d, MAXGEOMSPOINST %d", nelems, MAXGEOMSPOINTS);
-
- while (!result)
+ for (i=0; i<nelems; i++)
{
PG_LWGEOM *geom = (PG_LWGEOM *)(ARR_DATA_PTR(array)+offset);
offset += INTALIGN(VARSIZE(geom));
-
+ pgis_geom = geom;
/* Check is3d flag */
if ( TYPE_HASZ(geom->type) ) is3d = 1;
-
- /* Check SRID homogeneity */
- if ( ! i ) SRID = pglwgeom_getSRID(geom);
- else errorIfSRIDMismatch(SRID, pglwgeom_getSRID(geom));
-
- geoms[ngeoms] = g1 = POSTGIS2GEOS(geom);
-
- npoints += GEOSGetNumCoordinate(geoms[ngeoms]);
-
- ++ngeoms;
- ++i;
-
- POSTGIS_DEBUGF(4, "Loop %d, npoints: %d", i, npoints);
-
- /*
- * Maximum count of geometry points reached
- * or end of them, collect and buffer(0).
- */
- if ( (npoints>=MAXGEOMSPOINTS && ngeoms>1) || i==nelems)
+ /* Check SRID homogeneity and initialize geos result */
+ if ( ! i )
{
- POSTGIS_DEBUGF(4, " CHUNK (ngeoms:%d, npoints:%d, left:%d)",
- ngeoms, npoints, nelems-i);
-
- collection = GEOSMakeCollection(GEOS_GEOMETRYCOLLECTION,
- geoms, ngeoms);
-
- geos_result = GEOSBuffer(collection, 0, 0);
- if ( geos_result == NULL )
- {
- GEOSGeom_destroy(g1);
- lwerror("GEOS buffer() threw an error!");
- }
- GEOSGeom_destroy(collection);
-
- POSTGIS_DEBUG(4, " Buffer() executed");
-
- /*
- * If there are no other geoms in input
- * we've finished, otherwise we push
- * the result back on the input stack.
- */
- if ( i == nelems )
- {
- POSTGIS_DEBUGF(4, " Final result points: %d",
- GEOSGetNumCoordinate(geos_result));
-
- GEOSSetSRID(geos_result, SRID);
- result = GEOS2POSTGIS(geos_result, is3d);
- GEOSGeom_destroy(geos_result);
-
- POSTGIS_DEBUG(4, " Result computed");
-
- }
- else
- {
- geoms[0] = geos_result;
- ngeoms=1;
- npoints = GEOSGetNumCoordinate(geoms[0]);
-
- POSTGIS_DEBUGF(4, " Result pushed back on lwgeoms array (npoints:%d)", npoints);
- }
+ geos_result = POSTGIS2GEOS(geom);
+ SRID = pglwgeom_getSRID(geom);
+ continue;
}
+ else
+ {
+ errorIfSRIDMismatch(SRID, pglwgeom_getSRID(geom));
+ }
+ g1 = POSTGIS2GEOS(pgis_geom);
+ g2 = GEOSUnion(g1,geos_result);
+ if ( g2 == NULL )
+ {
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(geos_result);
+ elog(ERROR,"GEOS union() threw an error!");
+ }
+ GEOSGeom_destroy(g1);
+ GEOSGeom_destroy(geos_result);
+ geos_result = g2;
}
+ GEOSSetSRID(geos_result, SRID);
+ result = GEOS2POSTGIS(geos_result, is3d);
+ GEOSGeom_destroy(geos_result);
+ if ( result == NULL )
+ {
+ elog(ERROR, "GEOS2POSTGIS returned an error");
+ PG_RETURN_NULL(); /* never get here */
+ }
- /* compressType(result); */
-
PG_RETURN_POINTER(result);
-
}
-#endif /* def UNITE_USING_BUFFER */
-
/*
* select geomunion(
* 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
Modified: trunk/lwgeom/postgis.sql.in.c
===================================================================
--- trunk/lwgeom/postgis.sql.in.c 2009-01-21 21:55:13 UTC (rev 3546)
+++ trunk/lwgeom/postgis.sql.in.c 2009-01-21 23:19:50 UTC (rev 3547)
@@ -4054,6 +4054,12 @@
AS 'MODULE_PATHNAME','unite_garray'
LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable);
+CREATEFUNCTION _unite_garray_fast (geometry[])
+ RETURNS geometry
+ AS 'MODULE_PATHNAME','unite_garray_fast'
+ LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable);
+
+
-- Deprecation in 1.2.3
CREATE AGGREGATE GeomUnion (
sfunc = geom_accum,
@@ -4069,6 +4075,12 @@
stype = geometry[],
finalfunc = ST_unite_garray
);
+CREATE AGGREGATE ST_Union_Fast (
+ sfunc = ST_geom_accum,
+ basetype = geometry,
+ stype = geometry[],
+ finalfunc = _unite_garray_fast
+ );
-- Deprecation in 1.2.3
CREATEFUNCTION relate(geometry,geometry)
More information about the postgis-commits
mailing list