From 81b81d780aeeef97cbfea9a142d540595281d6d8 Mon Sep 17 00:00:00 2001
From: Sandro Santilli <strk@keybit.net>
Date: Fri, 27 Aug 2004 14:35:26 +0000
Subject: [PATCH] Added LWGEOM_EXPLODED structure definition and utility funx.
 Added collector() function and memcollect() aggregate. Still faulting...

git-svn-id: http://svn.osgeo.org/postgis/trunk@759 b70326c6-7e19-0410-871a-916f4a2858ee
---
 lwgeom/MISSING_OBJECTS          |   2 -
 lwgeom/TODO                     |   2 +-
 lwgeom/lwgeom.h                 |  36 ++-
 lwgeom/lwgeom_api.c             | 476 +++++++++++++++++++++++++++++++-
 lwgeom/lwgeom_functions_basic.c | 135 ++++++++-
 lwgeom/lwgeom_inout.c           |   2 +-
 lwgeom/lwpostgis.sql.in         |  11 +
 7 files changed, 656 insertions(+), 8 deletions(-)

diff --git a/lwgeom/MISSING_OBJECTS b/lwgeom/MISSING_OBJECTS
index 119bd0666..996013f3f 100644
--- a/lwgeom/MISSING_OBJECTS
+++ b/lwgeom/MISSING_OBJECTS
@@ -1,7 +1,6 @@
 # This is a list of objects still missing lwgeom support
 
 AGGREGATE: KEEPING AGGREGATE [collect(geometry)]
-AGGREGATE: KEEPING AGGREGATE [memcollect(geometry)]
 FNCAST: KEEPING FNCAST geometry(text) (see CAST)
 FUNC: KEEPING FUNCTION: [line_interpolate_point(geometry, double precision)]
 FUNC: KEEPING FUNCTION: [simplify(geometry, double precision)]
@@ -10,7 +9,6 @@ FUNC: KEEPING FUNCTION: [asbinary(geometry, text)]
 FUNC: KEEPING FUNCTION: [boundary(geometry)]
 FUNC: KEEPING FUNCTION: [box(geometry)]
 FUNC: KEEPING FUNCTION: [collect_garray(geometry[])]
-FUNC: KEEPING FUNCTION: [collector(geometry, geometry)]
 FUNC: KEEPING FUNCTION: [envelope(geometry)]
 FUNC: KEEPING FUNCTION: [equals(geometry, geometry)]
 FUNC: KEEPING FUNCTION: [expand(geometry, double precision)]
diff --git a/lwgeom/TODO b/lwgeom/TODO
index 3820ee28a..6d0a00abf 100644
--- a/lwgeom/TODO
+++ b/lwgeom/TODO
@@ -1,4 +1,4 @@
-- Decide what to do with bogus max_distance(geom, geom) function
+- Find the bug in memcollect/lwexploded_serialize/...
 
 - Check spheroid misure functions
 	o can't understand if 3d computation works since
diff --git a/lwgeom/lwgeom.h b/lwgeom/lwgeom.h
index 0b5870b77..0fd449474 100644
--- a/lwgeom/lwgeom.h
+++ b/lwgeom/lwgeom.h
@@ -1,5 +1,6 @@
 //lwgeom.h
 
+
 // basic API for handling the LWGEOM, BOX2DFLOAT4, LWPOINT, LWLINE, and LWPOLY.
 // See below for other support types like POINTARRAY and LWGEOM_INSPECTED
 
@@ -286,7 +287,7 @@ extern POINT2D lwpoint_getPoint2d(LWPOINT *point);
 extern POINT3D lwpoint_getPoint3d(LWPOINT *point);
 
 //find length of this serialized point
-extern uint32 lwpoint_findlength(char *serialized_line);
+extern uint32 lwpoint_findlength(char *serialized_point);
 
 //--------------------------------------------------------
 
@@ -353,6 +354,38 @@ typedef struct
 	char  **sub_geoms;    // list of pointers (into serialized_form) of the sub-geoms
 } LWGEOM_INSPECTED;
 
+/*
+ * This structure is intended to be used for geometry collection construction
+ */
+typedef struct
+{
+	int SRID;
+	int ndims;
+	int npoints;
+	char **points;
+	int nlines;
+	char **lines;
+	int npolys;
+	char **polys;
+} LWGEOM_EXPLODED;
+
+void pfree_exploded(LWGEOM_EXPLODED *exploded);
+
+// Returns a 'palloced' union of the two input exploded geoms.
+// Returns NULL if SRID or ndims do not match.
+LWGEOM_EXPLODED * lwexploded_sum(LWGEOM_EXPLODED *exp1, LWGEOM_EXPLODED *exp2);
+
+/*
+ * This function recursively scan the given serialized geometry
+ * and returns a list of _all_ subgeoms in it (deep-first)
+ */
+LWGEOM_EXPLODED *lwgeom_explode(char *serialized);
+
+// Serialize an LWGEOM_EXPLODED object.
+// SRID and ndims will be taken from exploded structure.
+// wantbbox will determine result bbox.
+char *lwexploded_serialize(LWGEOM_EXPLODED *exploded, int wantbbox);
+
 // note - for a simple type (ie. point), this will have sub_geom[0] = serialized_form.
 // for multi-geomtries sub_geom[0] will be a few bytes into the serialized form
 // This function just computes the length of each sub-object and pre-caches this info.
@@ -441,6 +474,7 @@ extern int lwgeom_seralizedformlength_inspected(LWGEOM_INSPECTED *inspected, int
 // get the SRID from the LWGEOM
 // none present => -1
 extern int lwgeom_getSRID(LWGEOM *lwgeom);
+extern int lwgeom_getsrid(char *serialized);
 extern LWGEOM *lwgeom_setSRID(LWGEOM *lwgeom, int32 newSRID);
 
 //get bounding box of LWGEOM (automatically calls the sub-geometries bbox generators)
diff --git a/lwgeom/lwgeom_api.c b/lwgeom/lwgeom_api.c
index 0d6e52a97..74924e5cd 100644
--- a/lwgeom/lwgeom_api.c
+++ b/lwgeom/lwgeom_api.c
@@ -1202,7 +1202,7 @@ uint32 lwpoint_findlength(char *serialized_point)
 
 	type = (unsigned char) serialized_point[0];
 
-	if ( lwgeom_getType(type) != POINTTYPE) return -9999;
+	if ( lwgeom_getType(type) != POINTTYPE) return 0;
 
 #ifdef DEBUG
 elog(NOTICE, "lwpoint_findlength called (%d)", result);
@@ -1291,7 +1291,10 @@ LWPOLY *lwpoly_deserialize(char *serialized_form)
 	int t;
 
 	if (serialized_form == NULL)
+	{
+		elog(ERROR, "lwpoly_deserialize called with NULL arg");
 		return NULL;
+	}
 
 	result = (LWPOLY*) palloc(sizeof(LWPOLY));
 
@@ -1301,7 +1304,11 @@ LWPOLY *lwpoly_deserialize(char *serialized_form)
 	loc = serialized_form;
 
 	if ( lwgeom_getType(type) != POLYGONTYPE)
+	{
+		elog(ERROR, "lwpoly_deserialize called with arg of type %d",
+			lwgeom_getType(type));
 		return NULL;
+	}
 
 
 	loc = serialized_form+1;
@@ -1671,6 +1678,8 @@ LWGEOM_INSPECTED *lwgeom_inspect(char *serialized_form)
 		result->ngeometries);
 #endif
 
+	if ( ! result->ngeometries ) return result;
+
 	sub_geoms = (char**) palloc(sizeof(char*) * result->ngeometries );
 	result->sub_geoms = sub_geoms;
 	sub_geoms[0] = loc;
@@ -2317,7 +2326,8 @@ BOX3D *lw_geom_getBB_inspected(LWGEOM_INSPECTED *inspected)
 
 void pfree_inspected(LWGEOM_INSPECTED *inspected)
 {
-	pfree(inspected->sub_geoms);
+	if ( inspected->ngeometries )
+		pfree(inspected->sub_geoms);
 	pfree(inspected);
 }
 
@@ -2482,6 +2492,23 @@ void printType(unsigned char type)
 	elog(NOTICE,"type 0x%x ==> hasBBOX=%i, hasSRID=%i, ndims=%i, type=%i",(unsigned int) type, lwgeom_hasBBOX(type), lwgeom_hasSRID(type),lwgeom_ndims(type), lwgeom_getType(type));
 }
 
+// get the SRID from the LWGEOM
+// none present => -1
+int lwgeom_getsrid(char *serialized)
+{
+	unsigned char type = serialized[0];
+	char *loc = serialized+1;
+
+	if ( ! lwgeom_hasSRID(type)) return -1;
+
+	if (lwgeom_hasBBOX(type))
+	{
+		loc += sizeof(BOX2DFLOAT4);
+	}
+
+	return get_int32(loc);
+}
+
 // get the SRID from the LWGEOM
 // none present => -1
 int lwgeom_getSRID(LWGEOM *lwgeom)
@@ -2612,3 +2639,448 @@ LWGEOM_construct(char *ser, int SRID, int wantbbox)
 
 	return result;
 }
+
+void
+pfree_exploded(LWGEOM_EXPLODED *exploded)
+{
+	if ( exploded->npoints )
+		pfree(exploded->points);
+	if ( exploded->nlines )
+		pfree(exploded->lines);
+	if ( exploded->npolys )
+		pfree(exploded->polys);
+	pfree(exploded);
+};
+
+/*
+ * This function recursively scan the given serialized geometry
+ * and returns a list of _all_ subgeoms in it (deep-first)
+ */
+LWGEOM_EXPLODED *
+lwgeom_explode(char *serialized)
+{
+	LWGEOM_INSPECTED *inspected = lwgeom_inspect(serialized);
+	LWGEOM_EXPLODED *subexploded, *result;
+	int i;
+
+	result = palloc(sizeof(LWGEOM_EXPLODED));
+	result->points = palloc(1);
+	result->lines = palloc(1);
+	result->polys = palloc(1);
+	result->npoints = 0;
+	result->nlines = 0;
+	result->npolys = 0;
+
+	if ( ! inspected->ngeometries )
+	{
+		pfree(result->points);
+		pfree(result->lines);
+		pfree(result->polys);
+		result->SRID = -1;
+		result->ndims = 0;
+		return result;
+	}
+
+	result->SRID = lwgeom_getsrid(serialized);
+	result->ndims = lwgeom_ndims(serialized[0]);
+
+	for (i=0; i<inspected->ngeometries; i++)
+	{
+
+		char *subgeom = inspected->sub_geoms[i];
+		int type = lwgeom_getType(subgeom[0]);
+
+		if ( type == POINTTYPE )
+		{
+			result->points = repalloc(result->points,
+				result->npoints+1);
+			result->points[result->npoints] = subgeom;
+			result->npoints++;
+			continue;
+		}
+
+		if ( type == LINETYPE )
+		{
+			result->lines = repalloc(result->lines,
+				result->nlines+1);
+			result->lines[result->nlines] = subgeom;
+			result->nlines++;
+			continue;
+		}
+
+		if ( type == POLYGONTYPE )
+		{
+			result->polys = repalloc(result->polys,
+				result->npolys+1);
+			result->polys[result->npolys] = subgeom;
+			result->npolys++;
+			continue;
+		}
+
+		// it's a multi geometry, recurse
+		subexploded = lwgeom_explode(subgeom);
+
+		// Re-allocate adding space for new exploded geoms
+		// (-1 because 1 was already allocated for the collection)
+		// Copy subgeom pointers from subexploded to current
+		// exploded.
+
+		if ( subexploded->npoints )
+		{
+			result->points = repalloc(result->points,
+				result->npoints+subexploded->npoints-1);
+
+			memcpy(result->points[result->npoints],
+				subexploded->points,
+				subexploded->npoints*sizeof(char *));
+
+			result->npoints += subexploded->npoints;
+		}
+
+		if ( subexploded->nlines )
+		{
+			result->lines = repalloc(result->lines,
+				result->nlines+subexploded->nlines-1);
+
+			memcpy(result->lines[result->nlines],
+				subexploded->lines,
+				subexploded->nlines*sizeof(char *));
+
+			result->nlines += subexploded->nlines;
+		}
+
+		if ( subexploded->npolys )
+		{
+			result->polys = repalloc(result->polys,
+				result->npolys+subexploded->npolys-1);
+
+			memcpy(result->polys[result->npolys],
+				subexploded->polys,
+				subexploded->npolys*sizeof(char *));
+
+			result->npolys += subexploded->npolys;
+		}
+
+		// release subexploded memory
+		pfree_exploded(subexploded);
+
+	}
+
+	return result;
+}
+
+// Returns a 'palloced' union of the two input exploded geoms
+// Returns NULL if SRID or ndims do not match.
+LWGEOM_EXPLODED *
+lwexploded_sum(LWGEOM_EXPLODED *exp1, LWGEOM_EXPLODED *exp2)
+{
+	LWGEOM_EXPLODED *expcoll;
+	char *loc;
+
+	if ( exp1->ndims != exp2->ndims ) return NULL;
+	if ( exp1->SRID != exp2->SRID ) return NULL;
+
+	expcoll = palloc(sizeof(LWGEOM_EXPLODED));
+
+	expcoll->npoints = exp1->npoints + exp2->npoints;
+	if ( expcoll->npoints ) {
+		expcoll->points = (char **)palloc(expcoll->npoints*sizeof(char *));
+		loc = (char *)&(expcoll->points[0]);
+		if ( exp1->npoints ) {
+			memcpy(loc, exp1->points,
+				exp1->npoints*sizeof(char *));
+			loc += exp1->npoints*sizeof(char *);
+		}
+		if ( exp2->npoints ) {
+			memcpy(loc, exp2->points,
+				exp2->npoints*sizeof(char *));
+		}
+	}
+
+	expcoll->nlines = exp1->nlines + exp2->nlines;
+	if ( expcoll->nlines ) {
+		expcoll->lines = palloc(expcoll->nlines*sizeof(char *));
+		loc = (char *)&(expcoll->lines[0]);
+		if ( exp1->nlines ) {
+			memcpy(loc, exp1->lines,
+				exp1->nlines*sizeof(char *));
+			loc += exp1->nlines*sizeof(char *);
+		}
+		if ( exp2->nlines ) {
+			memcpy(loc, exp2->lines,
+				exp2->nlines*sizeof(char *));
+		}
+	}
+
+	expcoll->npolys = exp1->npolys + exp2->npolys;
+	if ( expcoll->npolys ) {
+		expcoll->polys = palloc(expcoll->npolys*sizeof(char *));
+		loc = (char *)&(expcoll->polys[0]);
+		if ( exp1->npolys ) {
+			memcpy(loc, exp1->polys,
+				exp1->npolys*sizeof(char *));
+			loc += exp1->npolys*sizeof(char *);
+		}
+		if ( exp2->npolys ) {
+			memcpy(loc, exp2->polys,
+				exp2->npolys*sizeof(char *));
+		}
+	}
+
+	expcoll->ndims = exp1->ndims;
+	expcoll->SRID = exp1->SRID;
+
+	return expcoll;
+}
+
+
+/*
+ * Serialized a LWGEOM_EXPLODED structure
+ */
+char *
+lwexploded_serialize(LWGEOM_EXPLODED *exploded, int wantbbox)
+{
+	unsigned int size=0;
+	int i;
+	int ntypes = 0;
+	int ngeoms = 0;
+	char *result, *loc;
+	int outtype = 0;
+	LWPOLY *poly;
+	LWLINE *line;
+	LWPOINT *point;
+	BOX2DFLOAT4 *box2d;
+	BOX3D *box3d;
+	char *ser;
+
+	// find size of all geoms.
+	// If BBOX and SRID are included this size could be
+	// larger then needed, but that should not be a problem
+	for (i=0; i<exploded->npoints; i++)
+		size += lwpoint_findlength(exploded->points[i]);
+	for (i=0; i<exploded->nlines; i++)
+		size += lwline_findlength(exploded->lines[i]);
+	for (i=0; i<exploded->npolys; i++)
+		size += lwpoly_findlength(exploded->polys[i]);
+
+	if ( exploded->npoints )
+	{
+		ntypes++;
+		outtype = (exploded->npoints>1) ? MULTIPOINTTYPE : POINTTYPE;
+	}
+	if ( exploded->nlines )
+	{
+		ntypes++;
+		if ( outtype ) outtype = COLLECTIONTYPE;
+		else outtype = (exploded->nlines>1) ? MULTILINETYPE : LINETYPE;
+	}
+	if ( exploded->npolys )
+	{
+		ntypes++;
+		if ( outtype ) outtype = COLLECTIONTYPE;
+		else outtype = (exploded->npolys>1) ? MULTIPOLYGONTYPE : POLYGONTYPE;
+	}
+
+	ngeoms = exploded->npoints + exploded->nlines + exploded->npolys;
+
+#ifdef DEBUG
+	elog(NOTICE, " computed outtype: %d, ngeoms: %d", outtype, ngeoms);
+#endif
+
+
+	// For a single geometry just set SRID and BBOX (if requested)
+	if ( ngeoms < 2 )
+	{
+		if ( exploded->npoints )
+		{
+			point = lwpoint_deserialize(exploded->points[0]);
+			point->SRID = exploded->SRID;
+			ser = lwpoint_serialize(point);
+			pfree_point(point);
+			size = lwpoint_findlength(ser);
+		}
+		else if ( exploded->nlines )
+		{
+			line = lwline_deserialize(exploded->lines[0]);
+			line->SRID = exploded->SRID;
+			ser = lwline_serialize(line);
+			pfree_line(line);
+			size = lwline_findlength(ser);
+		}
+		else if ( exploded->npolys )
+		{
+			poly = lwpoly_deserialize(exploded->polys[0]);
+			poly->SRID = exploded->SRID;
+			ser = lwpoly_serialize(poly);
+			pfree_polygon(poly);
+			size = lwpoly_findlength(ser);
+		}
+		else return NULL;
+		if ( wantbbox && ! lwgeom_hasBBOX(ser[0]) )
+		{
+			result = palloc(size+4);
+			result[0] = TYPE_SETHASBBOX(ser[0], 1);
+			loc = result+1;
+			box3d = lw_geom_getBB_simple(ser);
+			box2d = box3d_to_box2df(box3d);
+			memcpy(loc, box2d, sizeof(BOX2DFLOAT4));
+			loc += sizeof(BOX2DFLOAT4);
+			memcpy(loc, (ser+1), size-1);
+			pfree(ser);
+			return result;
+		}
+		else
+		{
+			return ser;
+		}
+	}
+
+	// Add size for 3 multigeoms + root geom + bbox and srid.
+	// Also in this case we are considering worst case.
+	size += 24+sizeof(BOX2DFLOAT4); 
+
+#ifdef DEBUG
+	elog(NOTICE, " computed totsize: %d", size);
+#endif
+
+	result = palloc(size*2);
+	loc = result+1; // skip type
+
+	if ( wantbbox ) loc += sizeof(BOX2DFLOAT4); // skip box
+	if ( exploded->SRID != -1 ) loc += 4; // skip SRID
+
+	// If we have more then one type of geom
+	// write that number in the 'ngeoms' field of the
+	// output serialized form (internal geoms would be multi themself)
+	if ( ntypes > 1 )
+	{
+		memcpy(loc, &ntypes, 4);
+		loc += 4; 
+	}
+
+	else
+	{
+		loc--; // let the type be specified later.
+	}
+	
+	if ( exploded->npoints > 1 )
+	{
+		loc[0] = lwgeom_makeType_full(exploded->ndims, 0,
+			MULTIPOINTTYPE, 0);
+		loc++;
+		memcpy(loc, &exploded->npoints, 4); // numpoints
+		loc += 4; 
+	}
+	// Serialize points stripping BBOX and SRID if any
+	for (i=0; i<exploded->npoints; i++)
+	{
+		int subsize;
+
+		point = lwpoint_deserialize(exploded->points[i]);
+		point->SRID = -1;
+		ser = lwpoint_serialize(point);
+		subsize = lwpoint_findlength(ser);
+		memcpy(loc, ser, subsize);
+		loc += subsize;
+	}
+
+	if ( exploded->nlines > 1 )
+	{
+		loc[0] = lwgeom_makeType_full(exploded->ndims, 0,
+			MULTILINETYPE, 0);
+		loc++;
+		memcpy(loc, &exploded->nlines, 4); // numlines
+		loc += 4; 
+	}
+	// Serialize lines stripping BBOX and SRID if any
+	for (i=0; i<exploded->nlines; i++)
+	{
+		char *ser;
+		int subsize;
+
+		line = lwline_deserialize(exploded->lines[i]);
+		if ( line == NULL )
+		{
+	elog(ERROR, "Error deserializing %dnt line from exploded geom", i);
+	return NULL;
+		}
+		line->SRID = -1;
+		ser = lwline_serialize(line);
+		pfree_line(line);
+		subsize = lwline_findlength(ser);
+		memcpy(loc, ser, subsize);
+		pfree(ser);
+		loc += subsize;
+	}
+
+	if ( exploded->npolys > 1 )
+	{
+		loc[0] = lwgeom_makeType_full(exploded->ndims, 0,
+			MULTIPOLYGONTYPE, 0);
+		loc++;
+		memcpy(loc, &exploded->npolys, 4); // numpolys
+		loc += 4; 
+	}
+	// Serialize polys stripping BBOX and SRID if any
+	for (i=0; i<exploded->npolys; i++)
+	{
+		char *ser;
+		int subsize;
+
+		poly = lwpoly_deserialize(exploded->polys[i]);
+		if ( poly == NULL )
+		{
+	elog(ERROR, "Error deserializing %dnt polygon from exploded geom", i);
+	return NULL;
+		}
+		poly->SRID = -1;
+		ser = lwpoly_serialize(poly);
+		pfree_polygon(poly);
+		subsize = lwpoly_findlength(ser);
+#ifdef DEBUG
+		elog(NOTICE, "size of polygon %d: %d", i, subsize);
+#endif
+		memcpy(loc, ser, subsize);
+		pfree(ser);
+		loc += subsize;
+	}
+
+	// Ok. now we need to add type, SRID and bbox 
+	result[0] = lwgeom_makeType_full(exploded->ndims,
+		(exploded->SRID!=-1), outtype, wantbbox);
+	loc = result+1;
+
+	if ( wantbbox )
+	{
+		box3d = lw_geom_getBB_simple(result);
+		box2d = box3d_to_box2df(box3d);
+		memcpy(loc, box2d, sizeof(BOX2DFLOAT4));
+		loc += sizeof(BOX2DFLOAT4);
+	}
+
+	if ( exploded->SRID != -1 )
+	{
+		memcpy(loc, &(exploded->SRID), 4);
+		loc += 4; // useless.. we've finished
+	}
+
+#ifdef DEBUG
+	elog(NOTICE, "lwexploded_serialize finished");
+	elog(NOTICE, " type: %d", lwgeom_getType(result[0]));
+	elog(NOTICE, " SRID: %d", lwgeom_getsrid(result));
+	if ( lwgeom_hasBBOX(result[0]) )
+	{
+		{
+			BOX2DFLOAT4 boxbuf;
+			getbox2d_p(result, &boxbuf);
+			elog(NOTICE, " BBOX: %f,%f %f,%f",
+				boxbuf.xmin, boxbuf.ymin,
+				boxbuf.xmax, boxbuf.ymax);
+		}
+	}
+	elog(NOTICE, " numgeoms: %d", lwgeom_getnumgeometries(result));
+#endif
+		
+	return result;
+}
+
diff --git a/lwgeom/lwgeom_functions_basic.c b/lwgeom/lwgeom_functions_basic.c
index 0e614ebf2..feb836d8f 100644
--- a/lwgeom/lwgeom_functions_basic.c
+++ b/lwgeom/lwgeom_functions_basic.c
@@ -11,7 +11,7 @@
 
 #include "lwgeom.h"
 
-#define DEBUG
+//#define DEBUG
 
 Datum combine_box2d(PG_FUNCTION_ARGS);
 Datum LWGEOM_mem_size(PG_FUNCTION_ARGS);
@@ -34,11 +34,13 @@ Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS);
 Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS);
 Datum LWGEOM_translate(PG_FUNCTION_ARGS);
 Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS);
+Datum LWGEOM_collect(PG_FUNCTION_ARGS);
 
 // internal
 char * lwgeom_summary_recursive(char *serialized, int offset);
 int32 lwgeom_npoints_recursive(char *serialized);
 int32 lwgeom_nrings_recursive(char *serialized);
+void dump_lwexploded(LWGEOM_EXPLODED *exploded);
 
 // general utilities (might be moved in lwgeom_api.c)
 double lwgeom_polygon_area(LWPOLY *poly);
@@ -1892,3 +1894,134 @@ Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS)
 
 	PG_RETURN_BOOL(lwgeom_pt_inside_circle(pt, cx, cy, rr));
 }
+
+void
+dump_lwexploded(LWGEOM_EXPLODED *exploded)
+{
+	int i;
+
+	elog(NOTICE, "SRID=%d ndims=%d", exploded->SRID, exploded->ndims);
+	elog(NOTICE, "%d points, %d lines, %d polygons",
+		exploded->npoints, exploded->nlines, exploded->npolys);
+
+	for (i=0; i<exploded->npoints; i++)
+	{
+		elog(NOTICE, "Point%d @ %p", i, exploded->points[i]);
+	}
+
+	for (i=0; i<exploded->nlines; i++)
+	{
+		elog(NOTICE, "Line%d @ %p", i, exploded->lines[i]);
+	}
+
+	for (i=0; i<exploded->npolys; i++)
+	{
+		elog(NOTICE, "Poly%d @ %p", i, exploded->polys[i]);
+	}
+}
+
+// collect( geom, geom ) returns a geometry which contains
+// all the sub_objects from both of the argument geometries
+// returned geometry is the simplest possible, based on the types
+// of the colelct objects
+// ie. if all are of either X or multiX, then a multiX is returned.
+PG_FUNCTION_INFO_V1(LWGEOM_collect);
+Datum LWGEOM_collect(PG_FUNCTION_ARGS)
+{
+	Pointer geom1_ptr = PG_GETARG_POINTER(0);
+	Pointer geom2_ptr = PG_GETARG_POINTER(1);
+	LWGEOM *geom1, *geom2, *result;
+	LWGEOM_EXPLODED *exp1, *exp2, *expcoll;
+	char *serialized_result;
+	int wantbbox = 0;
+	int size;
+
+	// return null if both geoms are null
+	if ( (geom1_ptr == NULL) && (geom2_ptr == NULL) )
+	{
+		PG_RETURN_NULL();
+	}
+
+	// return a copy of the second geom if only first geom is null
+	if (geom1_ptr == NULL)
+	{
+		geom2 = (LWGEOM *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(1));
+		PG_RETURN_POINTER(geom2);
+	}
+
+	// return a copy of the first geom if only second geom is null
+	if (geom2_ptr == NULL)
+	{
+		geom1 = (LWGEOM *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
+		PG_RETURN_POINTER(geom1);
+	}
+
+	geom1 = (LWGEOM *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
+	geom2 = (LWGEOM *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(1));
+
+	if ( lwgeom_getSRID(geom1) != lwgeom_getSRID(geom2) )
+	{
+		elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
+		PG_RETURN_NULL();
+	}
+
+	exp1 = lwgeom_explode(SERIALIZED_FORM(geom1));
+	if ( exp1->npoints + exp1->nlines + exp1->npolys == 0 )
+	{
+		pfree(geom1);
+		pfree_exploded(exp1);
+		PG_RETURN_POINTER(geom2);
+	}
+
+	exp2 = lwgeom_explode(SERIALIZED_FORM(geom2));
+	if ( exp2->npoints + exp2->nlines + exp2->npolys == 0 )
+	{
+		pfree(geom2);
+		pfree_exploded(exp1);
+		pfree_exploded(exp2);
+		PG_RETURN_POINTER(geom1);
+	}
+
+	wantbbox = lwgeom_hasBBOX(geom1->type) ||
+		lwgeom_hasBBOX(geom2->type);
+
+	// Ok, make a new LWGEOM_EXPLODED being the union
+	// of the input ones
+
+	expcoll = lwexploded_sum(exp1, exp2);
+	if ( !expcoll ) 
+	{
+		elog(ERROR, "Could not sum exploded geoms");
+		PG_RETURN_NULL();
+	}
+
+	//DEBUG
+	dump_lwexploded(expcoll);
+
+	// Now serialized collected LWGEOM_EXPLODED
+	serialized_result = lwexploded_serialize(expcoll, wantbbox);
+	//serialized_result = lwexploded_serialize(exp1, wantbbox);
+	if ( ! serialized_result )
+	{
+		elog(ERROR, "Could not serialize exploded geoms");
+		PG_RETURN_NULL();
+	}
+
+	elog(NOTICE, "Serialized lwexploded"); 
+
+	// And create LWGEOM type (could provide a _buf version of
+	// the serializer instead)
+	size = lwgeom_seralizedformlength_simple(serialized_result);
+	result = LWGEOM_construct(serialized_result,
+		lwgeom_getsrid(serialized_result), wantbbox);
+	pfree(serialized_result);
+
+	pfree(geom1);
+	pfree(geom2);
+	pfree_exploded(exp1);
+	pfree_exploded(exp2);
+	pfree_exploded(expcoll);
+
+	//elog(ERROR, "Not implemented yet");
+	PG_RETURN_POINTER(result);
+}
diff --git a/lwgeom/lwgeom_inout.c b/lwgeom/lwgeom_inout.c
index 45ee1c9b4..ac542cba5 100644
--- a/lwgeom/lwgeom_inout.c
+++ b/lwgeom/lwgeom_inout.c
@@ -18,7 +18,7 @@
 #include "stringBuffer.h"
 
 
-#define DEBUG
+//#define DEBUG
 
 #include "lwgeom_pg.h"
 #include "wktparse.h"
diff --git a/lwgeom/lwpostgis.sql.in b/lwgeom/lwpostgis.sql.in
index 480583006..141e8a4a3 100644
--- a/lwgeom/lwpostgis.sql.in
+++ b/lwgeom/lwpostgis.sql.in
@@ -1186,6 +1186,17 @@ CREATEFUNCTION multi(geometry)
 	AS '@MODULE_FILENAME@', 'LWGEOM_force_multi'
 	LANGUAGE 'C' WITH (isstrict);
 
+CREATEFUNCTION collector(geometry, geometry) 
+	RETURNS geometry
+	AS '@MODULE_FILENAME@', 'LWGEOM_collect'
+	LANGUAGE 'C';
+
+CREATE AGGREGATE memcollect(
+	sfunc = collector,
+	basetype = geometry,
+	stype = geometry
+	);
+
 ------------------------------------------------------------------------
 
 --
-- 
2.40.0