Context Navigation

lwgeom_geos.c

Visit:

Revision 9817, 87.5 KB (checked in by strk, 2 days ago)
Do not assume geos allocates using malloc. Reduce memory use too.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision` Property svn:mime-type set to `text/plain`

Line
1	/**********************************************************************
2	*
3	* PostGIS - Spatial Types for PostgreSQL
4	* http://postgis.refractions.net
5	*
6	* Copyright 2009-2012 Sandro Santilli <strk@keybit.net>
7	* Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
8	* Copyright 2001-2003 Refractions Research Inc.
9	*
10	* This is free software; you can redistribute and/or modify it under
11	* the terms of the GNU General Public Licence. See the COPYING file.
12	*
13	**********************************************************************/
14
15	/* TODO: we probaby don't need _all_ these pgsql headers */
16	#include "postgres.h"
17	#include "fmgr.h"
18	#include "miscadmin.h"
19	#include "utils/array.h"
20	#include "utils/builtins.h"
21	#include "utils/hsearch.h"
22	#include "utils/memutils.h"
23	#include "executor/spi.h"
24	#include "funcapi.h"
25
26	#include "../postgis_config.h"
27	#include "lwgeom_functions_analytic.h" /* for point_in_polygon */
28	#include "lwgeom_cache.h"
29	#include "lwgeom_geos.h"
30	#include "liblwgeom_internal.h"
31	#include "lwgeom_rtree.h"
32
33	/*
34	** GEOS prepared geometry is only available from GEOS 3.1 onwards
35	*/
36	#define PREPARED_GEOM
37
38	#ifdef PREPARED_GEOM
39	#include "lwgeom_geos_prepared.h"
40	#endif
41
42	#include <string.h>
43	#include <assert.h>
44
45	/*
46	** Prototypes for SQL-bound functions
47	*/
48	Datum relate_full(PG_FUNCTION_ARGS);
49	Datum relate_pattern(PG_FUNCTION_ARGS);
50	Datum disjoint(PG_FUNCTION_ARGS);
51	Datum touches(PG_FUNCTION_ARGS);
52	Datum intersects(PG_FUNCTION_ARGS);
53	Datum crosses(PG_FUNCTION_ARGS);
54	Datum contains(PG_FUNCTION_ARGS);
55	Datum containsproperly(PG_FUNCTION_ARGS);
56	Datum covers(PG_FUNCTION_ARGS);
57	Datum overlaps(PG_FUNCTION_ARGS);
58	Datum isvalid(PG_FUNCTION_ARGS);
59	Datum isvalidreason(PG_FUNCTION_ARGS);
60	Datum isvaliddetail(PG_FUNCTION_ARGS);
61	Datum buffer(PG_FUNCTION_ARGS);
62	Datum intersection(PG_FUNCTION_ARGS);
63	Datum convexhull(PG_FUNCTION_ARGS);
64	Datum topologypreservesimplify(PG_FUNCTION_ARGS);
65	Datum difference(PG_FUNCTION_ARGS);
66	Datum boundary(PG_FUNCTION_ARGS);
67	Datum symdifference(PG_FUNCTION_ARGS);
68	Datum geomunion(PG_FUNCTION_ARGS);
69	Datum issimple(PG_FUNCTION_ARGS);
70	Datum isring(PG_FUNCTION_ARGS);
71	Datum pointonsurface(PG_FUNCTION_ARGS);
72	Datum GEOSnoop(PG_FUNCTION_ARGS);
73	Datum postgis_geos_version(PG_FUNCTION_ARGS);
74	Datum centroid(PG_FUNCTION_ARGS);
75	Datum polygonize_garray(PG_FUNCTION_ARGS);
76	Datum linemerge(PG_FUNCTION_ARGS);
77	Datum coveredby(PG_FUNCTION_ARGS);
78	Datum hausdorffdistance(PG_FUNCTION_ARGS);
79	Datum hausdorffdistancedensify(PG_FUNCTION_ARGS);
80	Datum ST_UnaryUnion(PG_FUNCTION_ARGS);
81	Datum ST_Equals(PG_FUNCTION_ARGS);
82	Datum ST_BuildArea(PG_FUNCTION_ARGS);
83
84	Datum pgis_union_geometry_array(PG_FUNCTION_ARGS);
85
86	/*
87	** Prototypes end
88	*/
89
90	static RTREE_POLY_CACHE *
91	GetRtreeCache(FunctionCallInfoData fcinfo, LWGEOM lwgeom, GSERIALIZED *poly)
92	{
93	MemoryContext old_context;
94	GeomCache* supercache = GetGeomCache(fcinfo);
95	RTREE_POLY_CACHE *poly_cache = supercache->rtree;
96
97	/*
98	* Switch the context to the function-scope context,
99	* retrieve the appropriate cache object, cache it for
100	* future use, then switch back to the local context.
101	*/
102	old_context = MemoryContextSwitchTo(fcinfo->flinfo->fn_mcxt);
103	poly_cache = retrieveCache(lwgeom, poly, poly_cache);
104	supercache->rtree = poly_cache;
105	MemoryContextSwitchTo(old_context);
106
107	return poly_cache;
108	}
109
110
111	PG_FUNCTION_INFO_V1(postgis_geos_version);
112	Datum postgis_geos_version(PG_FUNCTION_ARGS)
113	{
114	const char *ver = lwgeom_geos_version();
115	text *result = cstring2text(ver);
116	PG_RETURN_POINTER(result);
117	}
118
119
120	/**
121	* @brief Compute the Hausdorff distance thanks to the corresponding GEOS function
122	* @example hausdorffdistance {@link #hausdorffdistance} - SELECT st_hausdorffdistance(
123	* 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
124	* 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry);
125	*/
126
127	PG_FUNCTION_INFO_V1(hausdorffdistance);
128	Datum hausdorffdistance(PG_FUNCTION_ARGS)
129	{
130	#if POSTGIS_GEOS_VERSION < 32
131	lwerror("The GEOS version this PostGIS binary "
132	"was compiled against (%d) doesn't support "
133	"'ST_HausdorffDistance' function (3.2.0+ required)",
134	POSTGIS_GEOS_VERSION);
135	PG_RETURN_NULL();
136	#else
137	GSERIALIZED *geom1;
138	GSERIALIZED *geom2;
139	GEOSGeometry *g1;
140	GEOSGeometry *g2;
141	double result;
142	int retcode;
143
144	POSTGIS_DEBUG(2, "hausdorff_distance called");
145
146	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
147	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
148
149	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
150	PG_RETURN_NULL();
151
152	initGEOS(lwnotice, lwgeom_geos_error);
153
154	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
155	if ( 0 == g1 ) /* exception thrown at construction */
156	{
157	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
158	PG_RETURN_NULL();
159	}
160
161	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
162	if ( 0 == g2 ) /* exception thrown */
163	{
164	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
165	GEOSGeom_destroy(g1);
166	PG_RETURN_NULL();
167	}
168
169	retcode = GEOSHausdorffDistance(g1, g2, &result);
170	GEOSGeom_destroy(g1);
171	GEOSGeom_destroy(g2);
172
173	if (retcode == 0)
174	{
175	lwerror("GEOSHausdorffDistance: %s", lwgeom_geos_errmsg);
176	PG_RETURN_NULL(); /never get here /
177	}
178
179	PG_FREE_IF_COPY(geom1, 0);
180	PG_FREE_IF_COPY(geom2, 1);
181
182	PG_RETURN_FLOAT8(result);
183	#endif
184	}
185
186	/**
187	* @brief Compute the Hausdorff distance with densification thanks to the corresponding GEOS function
188	* @example hausdorffdistancedensify {@link #hausdorffdistancedensify} - SELECT st_hausdorffdistancedensify(
189	* 'POLYGON((0 0, 0 2, 1 2, 2 2, 2 0, 0 0))'::geometry,
190	* 'POLYGON((0.5 0.5, 0.5 2.5, 1.5 2.5, 2.5 2.5, 2.5 0.5, 0.5 0.5))'::geometry, 0.5);
191	*/
192
193	PG_FUNCTION_INFO_V1(hausdorffdistancedensify);
194	Datum hausdorffdistancedensify(PG_FUNCTION_ARGS)
195	{
196	#if POSTGIS_GEOS_VERSION < 32
197	lwerror("The GEOS version this PostGIS binary "
198	"was compiled against (%d) doesn't support "
199	"'ST_HausdorffDistance' function (3.2.0+ required)",
200	POSTGIS_GEOS_VERSION);
201	PG_RETURN_NULL();
202	#else
203	GSERIALIZED *geom1;
204	GSERIALIZED *geom2;
205	GEOSGeometry *g1;
206	GEOSGeometry *g2;
207	double densifyFrac;
208	double result;
209	int retcode;
210
211
212	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
213	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
214	densifyFrac = PG_GETARG_FLOAT8(2);
215
216	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
217	PG_RETURN_NULL();
218
219	initGEOS(lwnotice, lwgeom_geos_error);
220
221	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
222	if ( 0 == g1 ) /* exception thrown at construction */
223	{
224	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
225	PG_RETURN_NULL();
226	}
227
228	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
229	if ( 0 == g2 ) /* exception thrown at construction */
230	{
231	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
232	GEOSGeom_destroy(g1);
233	PG_RETURN_NULL();
234	}
235
236	retcode = GEOSHausdorffDistanceDensify(g1, g2, densifyFrac, &result);
237	GEOSGeom_destroy(g1);
238	GEOSGeom_destroy(g2);
239
240	if (retcode == 0)
241	{
242	lwerror("GEOSHausdorffDistanceDensify: %s", lwgeom_geos_errmsg);
243	PG_RETURN_NULL(); /never get here /
244	}
245
246	PG_FREE_IF_COPY(geom1, 0);
247	PG_FREE_IF_COPY(geom2, 1);
248
249	PG_RETURN_FLOAT8(result);
250	#endif
251	}
252
253
254
255	/**
256	* @brief This is the final function for GeomUnion
257	* aggregate. Will have as input an array of Geometries.
258	* Will iteratively call GEOSUnion on the GEOS-converted
259	* versions of them and return PGIS-converted version back.
260	* Changing combination order might speed up performance.
261	*/
262	PG_FUNCTION_INFO_V1(pgis_union_geometry_array);
263	Datum pgis_union_geometry_array(PG_FUNCTION_ARGS)
264	#if POSTGIS_GEOS_VERSION >= 33
265	{
266	/*
267	** For GEOS >= 3.3, use the new UnaryUnion functionality to merge the
268	** terminal collection from the ST_Union aggregate
269	*/
270	Datum datum;
271	ArrayType *array;
272
273	int is3d = LW_FALSE, gotsrid = LW_FALSE;
274	int nelems = 0, i = 0, geoms_size = 0, curgeom = 0;
275
276	GSERIALIZED *gser_out = NULL;
277
278	GEOSGeometry *g = NULL;
279	GEOSGeometry *g_union = NULL;
280	GEOSGeometry **geoms = NULL;
281
282	int srid = SRID_UNKNOWN;
283
284	size_t offset = 0;
285	bits8 *bitmap;
286	int bitmask;
287	int empty_type = 0;
288
289	datum = PG_GETARG_DATUM(0);
290
291	/* Null array, null geometry (should be empty?) */
292	if ( (Pointer *)datum == NULL ) PG_RETURN_NULL();
293
294	array = DatumGetArrayTypeP(datum);
295
296	/* How many things in our array? */
297	nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
298
299	/* PgSQL supplies a bitmap of which array entries are null */
300	bitmap = ARR_NULLBITMAP(array);
301
302	/* Empty array? Null return */
303	if ( nelems == 0 ) PG_RETURN_NULL();
304
305	/* One-element union is the element itself! */
306	if ( nelems == 1 )
307	{
308	/* If the element is a NULL then we need to handle it separately */
309	if (bitmap && (*bitmap & 1) == 0)
310	PG_RETURN_NULL();
311	else
312	PG_RETURN_POINTER((GSERIALIZED *)(ARR_DATA_PTR(array)));
313	}
314
315	/* Ok, we really need GEOS now ;) */
316	initGEOS(lwnotice, lwgeom_geos_error);
317
318	/*
319	** Collect the non-empty inputs and stuff them into a GEOS collection
320	*/
321	geoms_size = nelems;
322	geoms = palloc( sizeof(GEOSGeometry) geoms_size );
323
324	/*
325	** We need to convert the array of GSERIALIZED into a GEOS collection.
326	** First make an array of GEOS geometries.
327	*/
328	offset = 0;
329	bitmap = ARR_NULLBITMAP(array);
330	bitmask = 1;
331	for ( i = 0; i < nelems; i++ )
332	{
333	/* Only work on non-NULL entries in the array */
334	if ((bitmap && (*bitmap & bitmask) != 0) \|\| !bitmap)
335	{
336	GSERIALIZED gser_in = (GSERIALIZED )(ARR_DATA_PTR(array)+offset);
337
338	/* Check for SRID mismatch in array elements */
339	if ( gotsrid )
340	{
341	error_if_srid_mismatch(srid, gserialized_get_srid(gser_in));
342	}
343	else
344	{
345	/* Initialize SRID/dimensions info */
346	srid = gserialized_get_srid(gser_in);
347	is3d = gserialized_has_z(gser_in);
348	gotsrid = 1;
349	}
350
351	/* Don't include empties in the union */
352	if ( gserialized_is_empty(gser_in) )
353	{
354	int gser_type = gserialized_get_type(gser_in);
355	if (gser_type > empty_type)
356	{
357	empty_type = gser_type;
358	POSTGIS_DEBUGF(4, "empty_type = %d gser_type = %d", empty_type, gser_type);
359	}
360	}
361	else
362	{
363	g = (GEOSGeometry *)POSTGIS2GEOS(gser_in);
364
365	/* Uh oh! Exception thrown at construction... */
366	if ( ! g )
367	{
368	lwerror("One of the geometries in the set "
369	"could not be converted to GEOS: %s", lwgeom_geos_errmsg);
370	PG_RETURN_NULL();
371	}
372
373	/* Ensure we have enough space in our storage array */
374	if ( curgeom == geoms_size )
375	{
376	geoms_size *= 2;
377	geoms = repalloc( geoms, sizeof(GEOSGeometry) geoms_size );
378	}
379
380	geoms[curgeom] = g;
381	curgeom++;
382	}
383	offset += INTALIGN(VARSIZE(gser_in));
384	}
385
386	/* Advance NULL bitmap */
387	if (bitmap)
388	{
389	bitmask <<= 1;
390	if (bitmask == 0x100)
391	{
392	bitmap++;
393	bitmask = 1;
394	}
395	}
396	}
397
398	/*
399	** Take our GEOS geometries and turn them into a GEOS collection,
400	** then pass that into cascaded union.
401	*/
402	if (curgeom > 0)
403	{
404	g = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, geoms, curgeom);
405	if ( ! g )
406	{
407	lwerror("Could not create GEOS COLLECTION from geometry array: %s", lwgeom_geos_errmsg);
408	PG_RETURN_NULL();
409	}
410
411	g_union = GEOSUnaryUnion(g);
412	GEOSGeom_destroy(g);
413	if ( ! g_union )
414	{
415	lwerror("GEOSUnaryUnion: %s",
416	lwgeom_geos_errmsg);
417	PG_RETURN_NULL();
418	}
419
420	GEOSSetSRID(g_union, srid);
421	gser_out = GEOS2POSTGIS(g_union, is3d);
422	GEOSGeom_destroy(g_union);
423	}
424	/* No real geometries in our array, any empties? */
425	else
426	{
427	/* If it was only empties, we'll return the largest type number */
428	if ( empty_type > 0 )
429	{
430	PG_RETURN_POINTER(geometry_serialize(lwgeom_construct_empty(empty_type, srid, is3d, 0)));
431	}
432	/* Nothing but NULL, returns NULL */
433	else
434	{
435	PG_RETURN_NULL();
436	}
437	}
438
439	if ( ! gser_out )
440	{
441	/* Union returned a NULL geometry */
442	PG_RETURN_NULL();
443	}
444
445	PG_RETURN_POINTER(gser_out);
446	}
447
448	#else
449	{
450	/* For GEOS < 3.3, use the old CascadedUnion function for polygons and
451	brute force two-by-two for other types. */
452	Datum datum;
453	ArrayType *array;
454	int is3d = 0;
455	int nelems, i;
456	GSERIALIZED *result = NULL;
457	GSERIALIZED *pgis_geom = NULL;
458	GEOSGeometry * g1 = NULL;
459	GEOSGeometry * g2 = NULL;
460	GEOSGeometry * geos_result=NULL;
461	int srid=SRID_UNKNOWN;
462	size_t offset = 0;
463	bits8 *bitmap;
464	int bitmask;
465	#if POSTGIS_DEBUG_LEVEL > 0
466	static int call=1;
467	#endif
468	int gotsrid = 0;
469	int allpolys=1;
470
471	#if POSTGIS_DEBUG_LEVEL > 0
472	call++;
473	POSTGIS_DEBUGF(2, "GEOS incremental union (call %d)", call);
474	#endif
475
476	datum = PG_GETARG_DATUM(0);
477
478	/* TODO handle empties */
479
480	/* Null array, null geometry (should be empty?) */
481	if ( (Pointer *)datum == NULL ) PG_RETURN_NULL();
482
483	array = DatumGetArrayTypeP(datum);
484
485	nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
486
487	bitmap = ARR_NULLBITMAP(array);
488
489	POSTGIS_DEBUGF(3, "unite_garray: number of elements: %d", nelems);
490
491	if ( nelems == 0 ) PG_RETURN_NULL();
492
493	/* One-element union is the element itself */
494	if ( nelems == 1 )
495	{
496	/* If the element is a NULL then we need to handle it separately */
497	if (bitmap && (*bitmap & 1) == 0)
498	PG_RETURN_NULL();
499	else
500	PG_RETURN_POINTER((GSERIALIZED *)(ARR_DATA_PTR(array)));
501	}
502
503	/* Ok, we really need geos now ;) */
504	initGEOS(lwnotice, lwgeom_geos_error);
505
506	/*
507	** First, see if all our elements are POLYGON/MULTIPOLYGON
508	** If they are, we can use UnionCascaded for faster results.
509	*/
510	offset = 0;
511	bitmask = 1;
512	gotsrid = 0;
513	for ( i = 0; i < nelems; i++ )
514	{
515	/* Don't do anything for NULL values */
516	if ((bitmap && (*bitmap & bitmask) != 0) \|\| !bitmap)
517	{
518	GSERIALIZED pggeom = (GSERIALIZED )(ARR_DATA_PTR(array)+offset);
519	int pgtype = gserialized_get_type(pggeom);
520	offset += INTALIGN(VARSIZE(pggeom));
521	if ( ! gotsrid ) /* Initialize SRID */
522	{
523	srid = gserialized_get_srid(pggeom);
524	gotsrid = 1;
525	if ( gserialized_has_z(pggeom) ) is3d = 1;
526	}
527	else
528	{
529	error_if_srid_mismatch(srid, gserialized_get_srid(pggeom));
530	}
531
532	if ( pgtype != POLYGONTYPE && pgtype != MULTIPOLYGONTYPE )
533	{
534	allpolys = 0;
535	break;
536	}
537	}
538
539	/* Advance NULL bitmap */
540	if (bitmap)
541	{
542	bitmask <<= 1;
543	if (bitmask == 0x100)
544	{
545	bitmap++;
546	bitmask = 1;
547	}
548	}
549	}
550
551	if ( allpolys )
552	{
553	/*
554	** Everything is polygonal, so let's run the cascaded polygon union!
555	*/
556	int geoms_size = nelems;
557	int curgeom = 0;
558	GEOSGeometry **geoms = NULL;
559	geoms = palloc( sizeof(GEOSGeometry ) geoms_size );
560	/*
561	** We need to convert the array of GSERIALIZED into a GEOS MultiPolygon.
562	** First make an array of GEOS Polygons.
563	*/
564	offset = 0;
565	bitmap = ARR_NULLBITMAP(array);
566	bitmask = 1;
567	for ( i = 0; i < nelems; i++ )
568	{
569	GEOSGeometry* g;
570
571	/* Don't do anything for NULL values */
572	if ((bitmap && (*bitmap & bitmask) != 0) \|\| !bitmap)
573	{
574	GSERIALIZED pggeom = (GSERIALIZED )(ARR_DATA_PTR(array)+offset);
575	int pgtype = gserialized_get_type(pggeom);
576	offset += INTALIGN(VARSIZE(pggeom));
577	if ( pgtype == POLYGONTYPE )
578	{
579	if ( curgeom == geoms_size )
580	{
581	geoms_size *= 2;
582	geoms = repalloc( geoms, sizeof(GEOSGeom) * geoms_size );
583	}
584	g = (GEOSGeometry *)POSTGIS2GEOS(pggeom);
585	if ( 0 == g ) /* exception thrown at construction */
586	{
587	/* TODO: release GEOS allocated memory ! */
588	lwerror("One of the geometries in the set "
589	"could not be converted to GEOS: %s", lwgeom_geos_errmsg);
590	PG_RETURN_NULL();
591	}
592	geoms[curgeom] = g;
593	curgeom++;
594	}
595	if ( pgtype == MULTIPOLYGONTYPE )
596	{
597	int j = 0;
598	LWMPOLY lwmpoly = (LWMPOLY)lwgeom_from_gserialized(pggeom);;
599	for ( j = 0; j < lwmpoly->ngeoms; j++ )
600	{
601	GEOSGeometry* g;
602	if ( curgeom == geoms_size )
603	{
604	geoms_size *= 2;
605	geoms = repalloc( geoms, sizeof(GEOSGeom) * geoms_size );
606	}
607	/* This builds a LWPOLY on top of the serialized form */
608	g = LWGEOM2GEOS(lwpoly_as_lwgeom(lwmpoly->geoms[j]));
609	if ( 0 == g ) /* exception thrown at construction */
610	{
611	/* TODO: cleanup all GEOS memory */
612	lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
613	PG_RETURN_NULL();
614	}
615	geoms[curgeom++] = g;
616	}
617	lwmpoly_free(lwmpoly);
618	}
619	}
620
621	/* Advance NULL bitmap */
622	if (bitmap)
623	{
624	bitmask <<= 1;
625	if (bitmask == 0x100)
626	{
627	bitmap++;
628	bitmask = 1;
629	}
630	}
631	}
632	/*
633	** Take our GEOS Polygons and turn them into a GEOS MultiPolygon,
634	** then pass that into cascaded union.
635	*/
636	if (curgeom > 0)
637	{
638	g1 = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, curgeom);
639	if ( ! g1 )
640	{
641	/* TODO: cleanup geoms memory */
642	lwerror("Could not create MULTIPOLYGON from geometry array: %s", lwgeom_geos_errmsg);
643	PG_RETURN_NULL();
644	}
645	g2 = GEOSUnionCascaded(g1);
646	GEOSGeom_destroy(g1);
647	if ( ! g2 )
648	{
649	lwerror("GEOSUnionCascaded: %s",
650	lwgeom_geos_errmsg);
651	PG_RETURN_NULL();
652	}
653
654	GEOSSetSRID(g2, srid);
655	result = GEOS2POSTGIS(g2, is3d);
656	GEOSGeom_destroy(g2);
657	}
658	else
659	{
660	/* All we found were NULLs, so let's return NULL */
661	PG_RETURN_NULL();
662	}
663	}
664	else
665	{
666	/*
667	** Heterogeneous result, let's slog through this one union at a time.
668	*/
669	offset = 0;
670	bitmap = ARR_NULLBITMAP(array);
671	bitmask = 1;
672	for (i=0; i<nelems; i++)
673	{
674	/* Don't do anything for NULL values */
675	if ((bitmap && (*bitmap & bitmask) != 0) \|\| !bitmap)
676	{
677	GSERIALIZED geom = (GSERIALIZED )(ARR_DATA_PTR(array)+offset);
678	offset += INTALIGN(VARSIZE(geom));
679
680	pgis_geom = geom;
681
682	POSTGIS_DEBUGF(3, "geom %d @ %p", i, geom);
683
684	/* Check is3d flag */
685	if ( gserialized_has_z(geom) ) is3d = 1;
686
687	/* Check SRID homogeneity and initialize geos result */
688	if ( ! geos_result )
689	{
690	geos_result = (GEOSGeometry *)POSTGIS2GEOS(geom);
691	if ( 0 == geos_result ) /* exception thrown at construction */
692	{
693	lwerror("geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
694	PG_RETURN_NULL();
695	}
696	srid = gserialized_get_srid(geom);
697	POSTGIS_DEBUGF(3, "first geom is a %s", lwtype_name(gserialized_get_type(geom)));
698	}
699	else
700	{
701	error_if_srid_mismatch(srid, gserialized_get_srid(geom));
702
703	g1 = POSTGIS2GEOS(pgis_geom);
704	if ( 0 == g1 ) /* exception thrown at construction */
705	{
706	/* TODO: release GEOS allocated memory ! */
707	lwerror("First argument geometry could not be converted to GEOS: %s",
708	lwgeom_geos_errmsg);
709	PG_RETURN_NULL();
710	}
711
712	POSTGIS_DEBUGF(3, "unite_garray(%d): adding geom %d to union (%s)",
713	call, i, lwtype_name(gserialized_get_type(geom)));
714
715	g2 = GEOSUnion(g1, geos_result);
716	if ( g2 == NULL )
717	{
718	GEOSGeom_destroy((GEOSGeometry *)g1);
719	GEOSGeom_destroy((GEOSGeometry *)geos_result);
720	lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
721	}
722	GEOSGeom_destroy((GEOSGeometry *)g1);
723	GEOSGeom_destroy((GEOSGeometry *)geos_result);
724	geos_result = g2;
725	}
726	}
727
728	/* Advance NULL bitmap */
729	if (bitmap)
730	{
731	bitmask <<= 1;
732	if (bitmask == 0x100)
733	{
734	bitmap++;
735	bitmask = 1;
736	}
737	}
738
739	}
740
741	/* If geos_result is set then we found at least one non-NULL geometry */
742	if (geos_result)
743	{
744	GEOSSetSRID(geos_result, srid);
745	result = GEOS2POSTGIS(geos_result, is3d);
746	GEOSGeom_destroy(geos_result);
747	}
748	else
749	{
750	/* All we found were NULLs, so let's return NULL */
751	PG_RETURN_NULL();
752	}
753
754	}
755
756	if ( result == NULL )
757	{
758	/* Union returned a NULL geometry */
759	PG_RETURN_NULL();
760	}
761
762	PG_RETURN_POINTER(result);
763
764	}
765	#endif /* POSTGIS_GEOS_VERSION >= 33 */
766
767	/**
768	* @example ST_UnaryUnion {@link #geomunion} SELECT ST_UnaryUnion(
769	* 'POLYGON((0 0, 10 0, 0 10, 10 10, 0 0))'
770	* );
771	*
772	*/
773	PG_FUNCTION_INFO_V1(ST_UnaryUnion);
774	Datum ST_UnaryUnion(PG_FUNCTION_ARGS)
775	{
776	#if POSTGIS_GEOS_VERSION < 33
777	PG_RETURN_NULL();
778	lwerror("The GEOS version this PostGIS binary "
779	"was compiled against (%d) doesn't support "
780	"'GEOSUnaryUnion' function (3.3.0+ required)",
781	POSTGIS_GEOS_VERSION);
782	PG_RETURN_NULL();
783	#else /* POSTGIS_GEOS_VERSION >= 33 */
784	GSERIALIZED *geom1;
785	int is3d;
786	int srid;
787	GEOSGeometry g1, g3;
788	GSERIALIZED *result;
789
790	POSTGIS_DEBUG(2, "in ST_UnaryUnion");
791
792	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
793
794	/* UnaryUnion(empty) == (empty) */
795	if ( gserialized_is_empty(geom1) )
796	PG_RETURN_POINTER(geom1);
797
798	is3d = ( gserialized_has_z(geom1) );
799
800	srid = gserialized_get_srid(geom1);
801
802	initGEOS(lwnotice, lwgeom_geos_error);
803
804	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
805
806	if ( 0 == g1 ) /* exception thrown at construction */
807	{
808	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
809	PG_RETURN_NULL();
810	}
811
812	POSTGIS_DEBUGF(3, "g1=%s", GEOSGeomToWKT(g1));
813
814	g3 = GEOSUnaryUnion(g1);
815
816	POSTGIS_DEBUGF(3, "g3=%s", GEOSGeomToWKT(g3));
817
818	GEOSGeom_destroy(g1);
819
820	if (g3 == NULL)
821	{
822	lwerror("GEOSUnion: %s", lwgeom_geos_errmsg);
823	PG_RETURN_NULL(); /* never get here */
824	}
825
826
827	GEOSSetSRID(g3, srid);
828
829	result = GEOS2POSTGIS(g3, is3d);
830
831	GEOSGeom_destroy(g3);
832
833	if (result == NULL)
834	{
835	elog(ERROR, "ST_UnaryUnion failed converting GEOS result Geometry to PostGIS format");
836	PG_RETURN_NULL(); /never get here /
837	}
838
839	/* compressType(result); */
840
841	PG_FREE_IF_COPY(geom1, 0);
842
843	PG_RETURN_POINTER(result);
844	#endif /* POSTGIS_GEOS_VERSION >= 33 */
845	}
846
847
848	/**
849	* @example geomunion {@link #geomunion} SELECT geomunion(
850	* 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
851	* 'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))'
852	* );
853	*
854	*/
855	PG_FUNCTION_INFO_V1(geomunion);
856	Datum geomunion(PG_FUNCTION_ARGS)
857	{
858	GSERIALIZED *geom1;
859	GSERIALIZED *geom2;
860	GSERIALIZED *result;
861	LWGEOM lwgeom1, lwgeom2, *lwresult ;
862
863	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
864	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
865
866	lwgeom1 = lwgeom_from_gserialized(geom1) ;
867	lwgeom2 = lwgeom_from_gserialized(geom2) ;
868
869	lwresult = lwgeom_union(lwgeom1, lwgeom2) ;
870	result = geometry_serialize(lwresult) ;
871
872	lwgeom_free(lwgeom1) ;
873	lwgeom_free(lwgeom2) ;
874	lwgeom_free(lwresult) ;
875
876	PG_FREE_IF_COPY(geom1, 0);
877	PG_FREE_IF_COPY(geom2, 1);
878
879	PG_RETURN_POINTER(result);
880	}
881
882
883	/**
884	* @example symdifference {@link #symdifference} - SELECT symdifference(
885	* 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
886	* 'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
887	*/
888	PG_FUNCTION_INFO_V1(symdifference);
889	Datum symdifference(PG_FUNCTION_ARGS)
890	{
891	GSERIALIZED *geom1;
892	GSERIALIZED *geom2;
893	GSERIALIZED *result;
894	LWGEOM lwgeom1, lwgeom2, *lwresult ;
895
896	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
897	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
898
899	lwgeom1 = lwgeom_from_gserialized(geom1) ;
900	lwgeom2 = lwgeom_from_gserialized(geom2) ;
901
902	lwresult = lwgeom_symdifference(lwgeom1, lwgeom2) ;
903	result = geometry_serialize(lwresult) ;
904
905	lwgeom_free(lwgeom1) ;
906	lwgeom_free(lwgeom2) ;
907	lwgeom_free(lwresult) ;
908
909	PG_FREE_IF_COPY(geom1, 0);
910	PG_FREE_IF_COPY(geom2, 1);
911
912	PG_RETURN_POINTER(result);
913	}
914
915
916	PG_FUNCTION_INFO_V1(boundary);
917	Datum boundary(PG_FUNCTION_ARGS)
918	{
919	GSERIALIZED *geom1;
920	GEOSGeometry g1, g3;
921	GSERIALIZED *result;
922	int srid;
923
924
925	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
926
927	/* Empty.Boundary() == Empty */
928	if ( gserialized_is_empty(geom1) )
929	PG_RETURN_POINTER(geom1);
930
931	srid = gserialized_get_srid(geom1);
932
933	initGEOS(lwnotice, lwgeom_geos_error);
934
935	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1 );
936
937	if ( 0 == g1 ) /* exception thrown at construction */
938	{
939	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
940	PG_RETURN_NULL();
941	}
942
943	g3 = (GEOSGeometry *)GEOSBoundary(g1);
944
945	if (g3 == NULL)
946	{
947	GEOSGeom_destroy(g1);
948	lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg);
949	PG_RETURN_NULL(); /* never get here */
950	}
951
952	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
953
954	GEOSSetSRID(g3, srid);
955
956	result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
957
958	if (result == NULL)
959	{
960	GEOSGeom_destroy(g1);
961
962	GEOSGeom_destroy(g3);
963	elog(NOTICE,"GEOS2POSTGIS threw an error (result postgis geometry formation)!");
964	PG_RETURN_NULL(); /* never get here */
965	}
966
967	GEOSGeom_destroy(g1);
968	GEOSGeom_destroy(g3);
969
970	PG_FREE_IF_COPY(geom1, 0);
971
972	PG_RETURN_POINTER(result);
973	}
974
975	PG_FUNCTION_INFO_V1(convexhull);
976	Datum convexhull(PG_FUNCTION_ARGS)
977	{
978	GSERIALIZED *geom1;
979	GEOSGeometry g1, g3;
980	GSERIALIZED *result;
981	LWGEOM *lwout;
982	int srid;
983	GBOX bbox;
984
985	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
986
987	/* Empty.ConvexHull() == Empty */
988	if ( gserialized_is_empty(geom1) )
989	PG_RETURN_POINTER(geom1);
990
991	srid = gserialized_get_srid(geom1);
992
993	initGEOS(lwnotice, lwgeom_geos_error);
994
995	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
996
997	if ( 0 == g1 ) /* exception thrown at construction */
998	{
999	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1000	PG_RETURN_NULL();
1001	}
1002
1003	g3 = (GEOSGeometry *)GEOSConvexHull(g1);
1004	GEOSGeom_destroy(g1);
1005
1006	if (g3 == NULL)
1007	{
1008	lwerror("GEOSConvexHull: %s", lwgeom_geos_errmsg);
1009	PG_RETURN_NULL(); /* never get here */
1010	}
1011
1012	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1013
1014	GEOSSetSRID(g3, srid);
1015
1016	lwout = GEOS2LWGEOM(g3, gserialized_has_z(geom1));
1017	GEOSGeom_destroy(g3);
1018
1019	if (lwout == NULL)
1020	{
1021	elog(ERROR,"convexhull() failed to convert GEOS geometry to LWGEOM");
1022	PG_RETURN_NULL(); /* never get here */
1023	}
1024
1025	/* Copy input bbox if any */
1026	if ( gserialized_get_gbox_p(geom1, &bbox) )
1027	{
1028	/* Force the box to have the same dimensionality as the lwgeom */
1029	bbox.flags = lwout->flags;
1030	lwout->bbox = gbox_copy(&bbox);
1031	}
1032
1033	result = geometry_serialize(lwout);
1034	lwgeom_free(lwout);
1035
1036	if (result == NULL)
1037	{
1038	elog(ERROR,"GEOS convexhull() threw an error (result postgis geometry formation)!");
1039	PG_RETURN_NULL(); /* never get here */
1040	}
1041
1042	PG_FREE_IF_COPY(geom1, 0);
1043	PG_RETURN_POINTER(result);
1044	}
1045
1046	PG_FUNCTION_INFO_V1(topologypreservesimplify);
1047	Datum topologypreservesimplify(PG_FUNCTION_ARGS)
1048	{
1049	GSERIALIZED *geom1;
1050	double tolerance;
1051	GEOSGeometry g1, g3;
1052	GSERIALIZED *result;
1053
1054	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1055	tolerance = PG_GETARG_FLOAT8(1);
1056
1057	/* Empty.Simplify() == Empty */
1058	if ( gserialized_is_empty(geom1) )
1059	PG_RETURN_POINTER(geom1);
1060
1061	initGEOS(lwnotice, lwgeom_geos_error);
1062
1063	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
1064	if ( 0 == g1 ) /* exception thrown at construction */
1065	{
1066	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1067	PG_RETURN_NULL();
1068	}
1069
1070	g3 = GEOSTopologyPreserveSimplify(g1,tolerance);
1071	GEOSGeom_destroy(g1);
1072
1073	if (g3 == NULL)
1074	{
1075	lwerror("GEOSTopologyPreserveSimplify: %s", lwgeom_geos_errmsg);
1076	PG_RETURN_NULL(); /* never get here */
1077	}
1078
1079	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1080
1081	GEOSSetSRID(g3, gserialized_get_srid(geom1));
1082
1083	result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
1084	GEOSGeom_destroy(g3);
1085
1086	if (result == NULL)
1087	{
1088	elog(ERROR,"GEOS topologypreservesimplify() threw an error (result postgis geometry formation)!");
1089	PG_RETURN_NULL(); /* never get here */
1090	}
1091
1092	PG_FREE_IF_COPY(geom1, 0);
1093	PG_RETURN_POINTER(result);
1094	}
1095
1096	PG_FUNCTION_INFO_V1(buffer);
1097	Datum buffer(PG_FUNCTION_ARGS)
1098	{
1099	GSERIALIZED *geom1;
1100	double size;
1101	GEOSGeometry g1, g3;
1102	GSERIALIZED *result;
1103	int quadsegs = 8; /* the default */
1104	int nargs;
1105	enum
1106	{
1107	ENDCAP_ROUND = 1,
1108	ENDCAP_FLAT = 2,
1109	ENDCAP_SQUARE = 3
1110	};
1111	enum
1112	{
1113	JOIN_ROUND = 1,
1114	JOIN_MITRE = 2,
1115	JOIN_BEVEL = 3
1116	};
1117	static const double DEFAULT_MITRE_LIMIT = 5.0;
1118	static const int DEFAULT_ENDCAP_STYLE = ENDCAP_ROUND;
1119	static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1120
1121	double mitreLimit = DEFAULT_MITRE_LIMIT;
1122	int endCapStyle = DEFAULT_ENDCAP_STYLE;
1123	int joinStyle = DEFAULT_JOIN_STYLE;
1124	char *param;
1125	char *params = NULL;
1126	LWGEOM *lwg;
1127
1128	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1129	size = PG_GETARG_FLOAT8(1);
1130
1131	/* Empty.Buffer() == Empty[polygon] */
1132	if ( gserialized_is_empty(geom1) )
1133	{
1134	lwg = lwpoly_as_lwgeom(lwpoly_construct_empty(
1135	gserialized_get_srid(geom1),
1136	0, 0)); // buffer wouldn't give back z or m anyway
1137	PG_RETURN_POINTER(geometry_serialize(lwg));
1138	}
1139
1140	nargs = PG_NARGS();
1141
1142	initGEOS(lwnotice, lwgeom_geos_error);
1143
1144	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
1145	if ( 0 == g1 ) /* exception thrown at construction */
1146	{
1147	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1148	PG_RETURN_NULL();
1149	}
1150
1151	if (nargs > 2)
1152	{
1153	/* We strdup `cause we're going to modify it */
1154	params = pstrdup(PG_GETARG_CSTRING(2));
1155
1156	POSTGIS_DEBUGF(3, "Params: %s", params);
1157
1158	for (param=params; ; param=NULL)
1159	{
1160	char key, val;
1161	param = strtok(param, " ");
1162	if ( param == NULL ) break;
1163	POSTGIS_DEBUGF(3, "Param: %s", param);
1164
1165	key = param;
1166	val = strchr(key, '=');
1167	if ( val == NULL \|\| *(val+1) == '\0' )
1168	{
1169	lwerror("Missing value for buffer "
1170	"parameter %s", key);
1171	break;
1172	}
1173	*val = '\0';
1174	++val;
1175
1176	POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1177
1178	if ( !strcmp(key, "endcap") )
1179	{
1180	/* Supported end cap styles:
1181	* "round", "flat", "square"
1182	*/
1183	if ( !strcmp(val, "round") )
1184	{
1185	endCapStyle = ENDCAP_ROUND;
1186	}
1187	else if ( !strcmp(val, "flat") \|\|
1188	!strcmp(val, "butt") )
1189	{
1190	endCapStyle = ENDCAP_FLAT;
1191	}
1192	else if ( !strcmp(val, "square") )
1193	{
1194	endCapStyle = ENDCAP_SQUARE;
1195	}
1196	else
1197	{
1198	lwerror("Invalid buffer end cap "
1199	"style: %s (accept: "
1200	"'round', 'flat', 'butt' "
1201	"or 'square'"
1202	")", val);
1203	break;
1204	}
1205
1206	}
1207	else if ( !strcmp(key, "join") )
1208	{
1209	if ( !strcmp(val, "round") )
1210	{
1211	joinStyle = JOIN_ROUND;
1212	}
1213	else if ( !strcmp(val, "mitre") \|\|
1214	!strcmp(val, "miter") )
1215	{
1216	joinStyle = JOIN_MITRE;
1217	}
1218	else if ( !strcmp(val, "bevel") )
1219	{
1220	joinStyle = JOIN_BEVEL;
1221	}
1222	else
1223	{
1224	lwerror("Invalid buffer end cap "
1225	"style: %s (accept: "
1226	"'round', 'mitre', 'miter' "
1227	" or 'bevel'"
1228	")", val);
1229	break;
1230	}
1231	}
1232	else if ( !strcmp(key, "mitre_limit") \|\|
1233	!strcmp(key, "miter_limit") )
1234	{
1235	/* mitreLimit is a float */
1236	mitreLimit = atof(val);
1237	}
1238	else if ( !strcmp(key, "quad_segs") )
1239	{
1240	/* quadrant segments is an int */
1241	quadsegs = atoi(val);
1242	}
1243	else
1244	{
1245	lwerror("Invalid buffer parameter: %s (accept: "
1246	"'endcap', 'join', 'mitre_limit', "
1247	"'miter_limit and "
1248	"'quad_segs')", key);
1249	break;
1250	}
1251	}
1252
1253	pfree(params); /* was pstrduped */
1254
1255	POSTGIS_DEBUGF(3, "endCap:%d joinStyle:%d mitreLimit:%g",
1256	endCapStyle, joinStyle, mitreLimit);
1257
1258	}
1259
1260	#if POSTGIS_GEOS_VERSION >= 32
1261
1262	g3 = GEOSBufferWithStyle(g1, size, quadsegs, endCapStyle, joinStyle, mitreLimit);
1263	GEOSGeom_destroy(g1);
1264
1265	#else /* POSTGIS_GEOS_VERSION < 32 */
1266
1267	if ( mitreLimit != DEFAULT_MITRE_LIMIT \|\|
1268	endCapStyle != DEFAULT_ENDCAP_STYLE \|\|
1269	joinStyle != DEFAULT_JOIN_STYLE )
1270	{
1271	lwerror("The GEOS version this PostGIS binary "
1272	"was compiled against (%d) doesn't support "
1273	"specifying a mitre limit != %d or styles different "
1274	"from 'round' (needs 3.2 or higher)",
1275	DEFAULT_MITRE_LIMIT, POSTGIS_GEOS_VERSION);
1276	}
1277
1278	g3 = GEOSBuffer(g1,size,quadsegs);
1279	GEOSGeom_destroy(g1);
1280
1281	#endif /* POSTGIS_GEOS_VERSION < 32 */
1282
1283	if (g3 == NULL)
1284	{
1285	lwerror("GEOSBuffer: %s", lwgeom_geos_errmsg);
1286	PG_RETURN_NULL(); /* never get here */
1287	}
1288
1289	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
1290
1291	GEOSSetSRID(g3, gserialized_get_srid(geom1));
1292
1293	result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
1294	GEOSGeom_destroy(g3);
1295
1296	if (result == NULL)
1297	{
1298	elog(ERROR,"GEOS buffer() threw an error (result postgis geometry formation)!");
1299	PG_RETURN_NULL(); /* never get here */
1300	}
1301
1302	PG_FREE_IF_COPY(geom1, 0);
1303	PG_RETURN_POINTER(result);
1304	}
1305
1306	/*
1307	* Compute at offset curve to a line
1308	*/
1309	Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
1310	PG_FUNCTION_INFO_V1(ST_OffsetCurve);
1311	Datum ST_OffsetCurve(PG_FUNCTION_ARGS)
1312	{
1313	#if POSTGIS_GEOS_VERSION < 32
1314	lwerror("The GEOS version this PostGIS binary "
1315	"was compiled against (%d) doesn't support "
1316	"ST_OffsetCurve function "
1317	"(needs 3.2 or higher)",
1318	POSTGIS_GEOS_VERSION);
1319	PG_RETURN_NULL(); /* never get here */
1320	#else
1321
1322	GSERIALIZED *gser_input;
1323	GSERIALIZED *gser_result;
1324	LWGEOM *lwgeom_input;
1325	LWGEOM *lwgeom_result;
1326	double size;
1327	int quadsegs = 8; /* the default */
1328	int nargs;
1329
1330	enum
1331	{
1332	JOIN_ROUND = 1,
1333	JOIN_MITRE = 2,
1334	JOIN_BEVEL = 3
1335	};
1336
1337	static const double DEFAULT_MITRE_LIMIT = 5.0;
1338	static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1339	double mitreLimit = DEFAULT_MITRE_LIMIT;
1340	int joinStyle = DEFAULT_JOIN_STYLE;
1341	char *param = NULL;
1342	char *paramstr = NULL;
1343
1344	/* Read SQL arguments */
1345	nargs = PG_NARGS();
1346	gser_input = (GSERIALIZED*) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1347	size = PG_GETARG_FLOAT8(1);
1348
1349	/* Check for a useable type */
1350	if ( gserialized_get_type(gser_input) != LINETYPE )
1351	{
1352	lwerror("ST_OffsetCurve only works with LineStrings");
1353	PG_RETURN_NULL();
1354	}
1355
1356	/*
1357	* For distance == 0, just return the input.
1358	* Note that due to a bug, GEOS 3.3.0 would return EMPTY.
1359	* See http://trac.osgeo.org/geos/ticket/454
1360	*/
1361	if ( size == 0 )
1362	PG_RETURN_POINTER(gser_input);
1363
1364	/* Read the lwgeom, check for errors */
1365	lwgeom_input = lwgeom_from_gserialized(gser_input);
1366	if ( ! lwgeom_input )
1367	lwerror("ST_OffsetCurve: lwgeom_from_gserialized returned NULL");
1368
1369	/* For empty inputs, just echo them back */
1370	if ( lwgeom_is_empty(lwgeom_input) )
1371	PG_RETURN_POINTER(gser_input);
1372
1373	/* Process the optional arguments */
1374	if ( nargs > 2 )
1375	{
1376	text *wkttext = PG_GETARG_TEXT_P(2);
1377	paramstr = text2cstring(wkttext);
1378
1379	POSTGIS_DEBUGF(3, "paramstr: %s", paramstr);
1380
1381	for ( param=paramstr; ; param=NULL )
1382	{
1383	char key, val;
1384	param = strtok(param, " ");
1385	if ( param == NULL ) break;
1386	POSTGIS_DEBUGF(3, "Param: %s", param);
1387
1388	key = param;
1389	val = strchr(key, '=');
1390	if ( val == NULL \|\| *(val+1) == '\0' )
1391	{
1392	lwerror("ST_OffsetCurve: Missing value for buffer parameter %s", key);
1393	break;
1394	}
1395	*val = '\0';
1396	++val;
1397
1398	POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1399
1400	if ( !strcmp(key, "join") )
1401	{
1402	if ( !strcmp(val, "round") )
1403	{
1404	joinStyle = JOIN_ROUND;
1405	}
1406	else if ( !(strcmp(val, "mitre") && strcmp(val, "miter")) )
1407	{
1408	joinStyle = JOIN_MITRE;
1409	}
1410	else if ( ! strcmp(val, "bevel") )
1411	{
1412	joinStyle = JOIN_BEVEL;
1413	}
1414	else
1415	{
1416	lwerror("Invalid buffer end cap style: %s (accept: "
1417	"'round', 'mitre', 'miter' or 'bevel')", val);
1418	break;
1419	}
1420	}
1421	else if ( !strcmp(key, "mitre_limit") \|\|
1422	!strcmp(key, "miter_limit") )
1423	{
1424	/* mitreLimit is a float */
1425	mitreLimit = atof(val);
1426	}
1427	else if ( !strcmp(key, "quad_segs") )
1428	{
1429	/* quadrant segments is an int */
1430	quadsegs = atoi(val);
1431	}
1432	else
1433	{
1434	lwerror("Invalid buffer parameter: %s (accept: "
1435	"'join', 'mitre_limit', 'miter_limit and "
1436	"'quad_segs')", key);
1437	break;
1438	}
1439	}
1440	POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g", joinStyle, mitreLimit);
1441	pfree(paramstr); /* alloc'ed in text2cstring */
1442	}
1443
1444	lwgeom_result = lwgeom_offsetcurve(lwgeom_as_lwline(lwgeom_input), size, quadsegs, joinStyle, mitreLimit);
1445
1446	if (lwgeom_result == NULL)
1447	lwerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");
1448
1449	gser_result = gserialized_from_lwgeom(lwgeom_result, 0, 0);
1450	lwgeom_free(lwgeom_input);
1451	lwgeom_free(lwgeom_result);
1452	PG_RETURN_POINTER(gser_result);
1453
1454	#endif /* POSTGIS_GEOS_VERSION < 32 */
1455	}
1456
1457
1458	PG_FUNCTION_INFO_V1(intersection);
1459	Datum intersection(PG_FUNCTION_ARGS)
1460	{
1461	GSERIALIZED *geom1;
1462	GSERIALIZED *geom2;
1463	GSERIALIZED *result;
1464	LWGEOM lwgeom1, lwgeom2, *lwresult ;
1465
1466	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1467	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1468
1469	lwgeom1 = lwgeom_from_gserialized(geom1) ;
1470	lwgeom2 = lwgeom_from_gserialized(geom2) ;
1471
1472	lwresult = lwgeom_intersection(lwgeom1, lwgeom2) ;
1473	result = geometry_serialize(lwresult) ;
1474
1475	lwgeom_free(lwgeom1) ;
1476	lwgeom_free(lwgeom2) ;
1477	lwgeom_free(lwresult) ;
1478
1479	PG_FREE_IF_COPY(geom1, 0);
1480	PG_FREE_IF_COPY(geom2, 1);
1481
1482	PG_RETURN_POINTER(result);
1483	}
1484
1485	/**
1486	* @example difference {@link #difference} - SELECT difference(
1487	* 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))',
1488	* 'POLYGON((5 5, 15 5, 15 7, 5 7, 5 5))');
1489	*/
1490	PG_FUNCTION_INFO_V1(difference);
1491	Datum difference(PG_FUNCTION_ARGS)
1492	{
1493	GSERIALIZED *geom1;
1494	GSERIALIZED *geom2;
1495	GSERIALIZED *result;
1496	LWGEOM lwgeom1, lwgeom2, *lwresult ;
1497
1498	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1499	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1500
1501	lwgeom1 = lwgeom_from_gserialized(geom1) ;
1502	lwgeom2 = lwgeom_from_gserialized(geom2) ;
1503
1504	lwresult = lwgeom_difference(lwgeom1, lwgeom2) ;
1505	result = geometry_serialize(lwresult) ;
1506
1507	lwgeom_free(lwgeom1) ;
1508	lwgeom_free(lwgeom2) ;
1509	lwgeom_free(lwresult) ;
1510
1511	PG_FREE_IF_COPY(geom1, 0);
1512	PG_FREE_IF_COPY(geom2, 1);
1513
1514	PG_RETURN_POINTER(result);
1515	}
1516
1517
1518	/**
1519	@example pointonsurface - {@link #pointonsurface} SELECT pointonsurface('POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))');
1520	*/
1521	PG_FUNCTION_INFO_V1(pointonsurface);
1522	Datum pointonsurface(PG_FUNCTION_ARGS)
1523	{
1524	GSERIALIZED *geom1;
1525	GEOSGeometry g1, g3;
1526	GSERIALIZED *result;
1527
1528	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1529
1530	/* Empty.PointOnSurface == Empty */
1531	if ( gserialized_is_empty(geom1) )
1532	PG_RETURN_POINTER(geom1);
1533
1534	initGEOS(lwnotice, lwgeom_geos_error);
1535
1536	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
1537
1538	if ( 0 == g1 ) /* exception thrown at construction */
1539	{
1540	elog(WARNING, "GEOSPointOnSurface(): %s", lwgeom_geos_errmsg);
1541	PG_RETURN_NULL();
1542	}
1543
1544	g3 = GEOSPointOnSurface(g1);
1545
1546	if (g3 == NULL)
1547	{
1548	GEOSGeom_destroy(g1);
1549	lwerror("GEOSPointOnSurface: %s", lwgeom_geos_errmsg);
1550	PG_RETURN_NULL(); /* never get here */
1551	}
1552
1553	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
1554
1555	GEOSSetSRID(g3, gserialized_get_srid(geom1));
1556
1557	result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
1558
1559	if (result == NULL)
1560	{
1561	GEOSGeom_destroy(g1);
1562	GEOSGeom_destroy(g3);
1563	elog(ERROR,"GEOS pointonsurface() threw an error (result postgis geometry formation)!");
1564	PG_RETURN_NULL(); /* never get here */
1565	}
1566
1567	GEOSGeom_destroy(g1);
1568	GEOSGeom_destroy(g3);
1569
1570	PG_FREE_IF_COPY(geom1, 0);
1571
1572	PG_RETURN_POINTER(result);
1573	}
1574
1575	PG_FUNCTION_INFO_V1(centroid);
1576	Datum centroid(PG_FUNCTION_ARGS)
1577	{
1578	GSERIALIZED geom, result;
1579	GEOSGeometry geosgeom, geosresult;
1580
1581	geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1582
1583	/* Empty.Centroid() == Empty */
1584	if ( gserialized_is_empty(geom) )
1585	PG_RETURN_POINTER(geom);
1586
1587	initGEOS(lwnotice, lwgeom_geos_error);
1588
1589	geosgeom = (GEOSGeometry *)POSTGIS2GEOS(geom);
1590
1591	if ( 0 == geosgeom ) /* exception thrown at construction */
1592	{
1593	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1594	PG_RETURN_NULL();
1595	}
1596
1597	geosresult = GEOSGetCentroid(geosgeom);
1598
1599	if ( geosresult == NULL )
1600	{
1601	GEOSGeom_destroy(geosgeom);
1602	lwerror("GEOSGetCentroid: %s", lwgeom_geos_errmsg);
1603	PG_RETURN_NULL();
1604	}
1605
1606	GEOSSetSRID(geosresult, gserialized_get_srid(geom));
1607
1608	result = GEOS2POSTGIS(geosresult, gserialized_has_z(geom));
1609
1610	if (result == NULL)
1611	{
1612	GEOSGeom_destroy(geosgeom);
1613	GEOSGeom_destroy(geosresult);
1614	elog(ERROR,"Error in GEOS-PGIS conversion");
1615	PG_RETURN_NULL();
1616	}
1617	GEOSGeom_destroy(geosgeom);
1618	GEOSGeom_destroy(geosresult);
1619
1620	PG_FREE_IF_COPY(geom, 0);
1621
1622	PG_RETURN_POINTER(result);
1623	}
1624
1625
1626
1627	/---------------------------------------------/
1628
1629
1630	/**
1631	* @brief Throws an ereport ERROR if either geometry is a COLLECTIONTYPE. Additionally
1632	* displays a HINT of the first 80 characters of the WKT representation of the
1633	* problematic geometry so a user knows which parameter and which geometry
1634	* is causing the problem.
1635	*/
1636	void errorIfGeometryCollection(GSERIALIZED g1, GSERIALIZED g2)
1637	{
1638	int t1 = gserialized_get_type(g1);
1639	int t2 = gserialized_get_type(g2);
1640
1641	char *hintmsg;
1642	char *hintwkt;
1643	size_t hintsz;
1644	LWGEOM *lwgeom;
1645
1646	if ( t1 == COLLECTIONTYPE)
1647	{
1648	lwgeom = lwgeom_from_gserialized(g1);
1649	hintwkt = lwgeom_to_wkt(lwgeom, WKT_SFSQL, DBL_DIG, &hintsz);
1650	hintmsg = lwmessage_truncate(hintwkt, 0, hintsz-1, 80, 1);
1651	ereport(ERROR,
1652	(errmsg("Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported."),
1653	errhint("Change argument 1: '%s'", hintmsg))
1654	);
1655	pfree(hintwkt);
1656	pfree(hintmsg);
1657	lwgeom_free(lwgeom);
1658	}
1659	else if (t2 == COLLECTIONTYPE)
1660	{
1661	lwgeom = lwgeom_from_gserialized(g2);
1662	hintwkt = lwgeom_to_wkt(lwgeom, WKT_SFSQL, DBL_DIG, &hintsz);
1663	hintmsg = lwmessage_truncate(hintwkt, 0, hintsz-1, 80, 1);
1664	ereport(ERROR,
1665	(errmsg("Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported."),
1666	errhint("Change argument 2: '%s'", hintmsg))
1667	);
1668	pfree(hintwkt);
1669	pfree(hintmsg);
1670	lwgeom_free(lwgeom);
1671	}
1672	}
1673
1674	PG_FUNCTION_INFO_V1(isvalid);
1675	Datum isvalid(PG_FUNCTION_ARGS)
1676	{
1677	GSERIALIZED *geom1;
1678	LWGEOM *lwgeom;
1679	bool result;
1680	GEOSGeom g1;
1681	#if POSTGIS_GEOS_VERSION < 33
1682	GBOX box1;
1683	#endif
1684
1685	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1686
1687	/* Empty.IsValid() == TRUE */
1688	if ( gserialized_is_empty(geom1) )
1689	PG_RETURN_BOOL(true);
1690
1691	#if POSTGIS_GEOS_VERSION < 33
1692	/* Short circuit and return FALSE if we have infinite coordinates */
1693	/* GEOS 3.3+ is supposed to handle this stuff OK */
1694	if ( gserialized_get_gbox_p(geom1, &box1) )
1695	{
1696	if ( isinf(box1.xmax) \|\| isinf(box1.ymax) \|\| isinf(box1.xmin) \|\| isinf(box1.ymin) \|\|
1697	isnan(box1.xmax) \|\| isnan(box1.ymax) \|\| isnan(box1.xmin) \|\| isnan(box1.ymin) )
1698	{
1699	lwnotice("Geometry contains an Inf or NaN coordinate");
1700	PG_RETURN_BOOL(FALSE);
1701	}
1702	}
1703	#endif
1704
1705	initGEOS(lwnotice, lwgeom_geos_error);
1706
1707	lwgeom = lwgeom_from_gserialized(geom1);
1708	if ( ! lwgeom )
1709	{
1710	lwerror("unable to deserialize input");
1711	}
1712	g1 = LWGEOM2GEOS(lwgeom);
1713	lwgeom_free(lwgeom);
1714
1715	if ( ! g1 )
1716	{
1717	/* should we drop the following
1718	* notice now that we have ST_isValidReason ?
1719	*/
1720	lwnotice("%s", lwgeom_geos_errmsg);
1721	PG_RETURN_BOOL(FALSE);
1722	}
1723
1724	result = GEOSisValid(g1);
1725	GEOSGeom_destroy(g1);
1726
1727	if (result == 2)
1728	{
1729	elog(ERROR,"GEOS isvalid() threw an error!");
1730	PG_RETURN_NULL(); /never get here /
1731	}
1732
1733	PG_FREE_IF_COPY(geom1, 0);
1734	PG_RETURN_BOOL(result);
1735	}
1736
1737	/*
1738	** IsValidReason is only available in the GEOS
1739	** C API > version 3.0
1740	*/
1741	PG_FUNCTION_INFO_V1(isvalidreason);
1742	Datum isvalidreason(PG_FUNCTION_ARGS)
1743	{
1744	GSERIALIZED *geom = NULL;
1745	char *reason_str = NULL;
1746	text *result = NULL;
1747	const GEOSGeometry *g1 = NULL;
1748	#if POSTGIS_GEOS_VERSION < 33
1749	GBOX box;
1750	#endif
1751
1752	geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1753
1754	#if POSTGIS_GEOS_VERSION < 33
1755	/* Short circuit and return if we have infinite coordinates */
1756	/* GEOS 3.3+ is supposed to handle this stuff OK */
1757	if ( gserialized_get_gbox_p(geom, &box) )
1758	{
1759	if ( isinf(box.xmax) \|\| isinf(box.ymax) \|\| isinf(box.xmin) \|\| isinf(box.ymin) \|\|
1760	isnan(box.xmax) \|\| isnan(box.ymax) \|\| isnan(box.xmin) \|\| isnan(box.ymin) )
1761	{
1762	const char *rsn = "Geometry contains an Inf or NaN coordinate";
1763	size_t len = strlen(rsn);
1764	result = palloc(VARHDRSZ + len);
1765	SET_VARSIZE(result, VARHDRSZ + len);
1766	memcpy(VARDATA(result), rsn, len);
1767	PG_FREE_IF_COPY(geom, 0);
1768	PG_RETURN_POINTER(result);
1769	}
1770	}
1771	#endif
1772
1773	initGEOS(lwnotice, lwgeom_geos_error);
1774
1775	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom);
1776	if ( g1 )
1777	{
1778	reason_str = GEOSisValidReason(g1);
1779	GEOSGeom_destroy((GEOSGeometry *)g1);
1780	if (reason_str == NULL)
1781	{
1782	elog(ERROR,"GEOSisValidReason() threw an error: %s", lwgeom_geos_errmsg);
1783	PG_RETURN_NULL(); /* never get here */
1784	}
1785	result = cstring2text(reason_str);
1786	GEOSFree(reason_str);
1787	}
1788	else
1789	{
1790	result = cstring2text(lwgeom_geos_errmsg);
1791	}
1792
1793
1794	PG_FREE_IF_COPY(geom, 0);
1795	PG_RETURN_POINTER(result);
1796
1797	}
1798
1799	/*
1800	** IsValidDetail is only available in the GEOS
1801	** C API >= version 3.3
1802	*/
1803	PG_FUNCTION_INFO_V1(isvaliddetail);
1804	Datum isvaliddetail(PG_FUNCTION_ARGS)
1805	{
1806	#if POSTGIS_GEOS_VERSION < 33
1807	lwerror("The GEOS version this PostGIS binary "
1808	"was compiled against (%d) doesn't support "
1809	"'isValidDetail' function (3.3.0+ required)",
1810	POSTGIS_GEOS_VERSION);
1811	PG_RETURN_NULL();
1812	#else /* POSTGIS_GEOS_VERSION >= 33 */
1813
1814	GSERIALIZED *geom = NULL;
1815	const GEOSGeometry *g1 = NULL;
1816	char values[3]; / valid bool, reason text, location geometry */
1817	char *geos_reason = NULL;
1818	char *reason = NULL;
1819	GEOSGeometry *geos_location = NULL;
1820	LWGEOM *location = NULL;
1821	char valid = 0;
1822	Datum result;
1823	TupleDesc tupdesc;
1824	HeapTuple tuple;
1825	AttInMetadata *attinmeta;
1826	int flags = 0;
1827
1828	/*
1829	* Build a tuple description for a
1830	* valid_detail tuple
1831	*/
1832	tupdesc = RelationNameGetTupleDesc("valid_detail");
1833	if ( ! tupdesc )
1834	{
1835	lwerror("TYPE valid_detail not found");
1836	PG_RETURN_NULL();
1837	}
1838
1839	/*
1840	* generate attribute metadata needed later to produce
1841	* tuples from raw C strings
1842	*/
1843	attinmeta = TupleDescGetAttInMetadata(tupdesc);
1844
1845	geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1846
1847	if ( PG_NARGS() > 1 && ! PG_ARGISNULL(1) ) {
1848	flags = PG_GETARG_INT32(1);
1849	}
1850
1851	initGEOS(lwnotice, lwgeom_geos_error);
1852
1853	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom);
1854
1855	if ( g1 )
1856	{
1857	valid = GEOSisValidDetail(g1, flags,
1858	&geos_reason, &geos_location);
1859	GEOSGeom_destroy((GEOSGeometry *)g1);
1860	if ( geos_reason )
1861	{
1862	reason = pstrdup(geos_reason);
1863	GEOSFree(geos_reason);
1864	}
1865	if ( geos_location )
1866	{
1867	location = GEOS2LWGEOM(geos_location, GEOSHasZ(geos_location));
1868	GEOSGeom_destroy((GEOSGeometry *)geos_location);
1869	}
1870
1871	if (valid == 2)
1872	{
1873	/* NOTE: should only happen on OOM or similar */
1874	lwerror("GEOS isvaliddetail() threw an exception!");
1875	PG_RETURN_NULL(); /* never gets here */
1876	}
1877	}
1878	else
1879	{
1880	/* TODO: check lwgeom_geos_errmsg for validity error */
1881	reason = pstrdup(lwgeom_geos_errmsg);
1882	}
1883
1884	/* the boolean validity */
1885	values[0] = valid ? "t" : "f";
1886
1887	/* the reason */
1888	values[1] = reason;
1889
1890	/* the location */
1891	values[2] = location ? lwgeom_to_hexwkb(location, WKB_EXTENDED, 0) : 0;
1892
1893	tuple = BuildTupleFromCStrings(attinmeta, values);
1894	result = HeapTupleGetDatum(tuple);
1895
1896	PG_RETURN_HEAPTUPLEHEADER(result);
1897
1898	#endif /* POSTGIS_GEOS_VERSION >= 33 */
1899	}
1900
1901	/**
1902	* overlaps(GSERIALIZED g1,GSERIALIZED g2)
1903	* @param g1
1904	* @param g2
1905	* @return if GEOS::g1->overlaps(g2) returns true
1906	* @throw an error (elog(ERROR,...)) if GEOS throws an error
1907	*/
1908	PG_FUNCTION_INFO_V1(overlaps);
1909	Datum overlaps(PG_FUNCTION_ARGS)
1910	{
1911	GSERIALIZED *geom1;
1912	GSERIALIZED *geom2;
1913	GEOSGeometry g1, g2;
1914	bool result;
1915	GBOX box1, box2;
1916
1917	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1918	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1919
1920	errorIfGeometryCollection(geom1,geom2);
1921	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
1922
1923	/* A.Overlaps(Empty) == FALSE */
1924	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
1925	PG_RETURN_BOOL(false);
1926
1927	/*
1928	* short-circuit 1: if geom2 bounding box does not overlap
1929	* geom1 bounding box we can prematurely return FALSE.
1930	* Do the test IFF BOUNDING BOX AVAILABLE.
1931	*/
1932	if ( gserialized_get_gbox_p(geom1, &box1) &&
1933	gserialized_get_gbox_p(geom2, &box2) )
1934	{
1935	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
1936	{
1937	PG_RETURN_BOOL(FALSE);
1938	}
1939	}
1940
1941	initGEOS(lwnotice, lwgeom_geos_error);
1942
1943	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
1944	if ( 0 == g1 ) /* exception thrown at construction */
1945	{
1946	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1947	PG_RETURN_NULL();
1948	}
1949
1950	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
1951
1952	if ( 0 == g2 ) /* exception thrown at construction */
1953	{
1954	GEOSGeom_destroy(g1);
1955	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
1956	PG_RETURN_NULL();
1957	}
1958
1959	result = GEOSOverlaps(g1,g2);
1960
1961	GEOSGeom_destroy(g1);
1962	GEOSGeom_destroy(g2);
1963	if (result == 2)
1964	{
1965	lwerror("GEOSOverlaps: %s", lwgeom_geos_errmsg);
1966	PG_RETURN_NULL(); /* never get here */
1967	}
1968
1969	PG_FREE_IF_COPY(geom1, 0);
1970	PG_FREE_IF_COPY(geom2, 1);
1971
1972	PG_RETURN_BOOL(result);
1973	}
1974
1975
1976	PG_FUNCTION_INFO_V1(contains);
1977	Datum contains(PG_FUNCTION_ARGS)
1978	{
1979	GSERIALIZED *geom1;
1980	GSERIALIZED *geom2;
1981	GEOSGeometry g1, g2;
1982	GBOX box1, box2;
1983	int type1, type2;
1984	LWGEOM *lwgeom;
1985	LWPOINT *point;
1986	RTREE_POLY_CACHE *poly_cache;
1987	bool result;
1988	#ifdef PREPARED_GEOM
1989	PrepGeomCache *prep_cache;
1990	#endif
1991
1992	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1993	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1994
1995	errorIfGeometryCollection(geom1,geom2);
1996	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
1997
1998	/* A.Contains(Empty) == FALSE */
1999	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2000	PG_RETURN_BOOL(false);
2001
2002	POSTGIS_DEBUG(3, "contains called.");
2003
2004	/*
2005	** short-circuit 1: if geom2 bounding box is not completely inside
2006	** geom1 bounding box we can prematurely return FALSE.
2007	** Do the test IFF BOUNDING BOX AVAILABLE.
2008	*/
2009	if ( gserialized_get_gbox_p(geom1, &box1) &&
2010	gserialized_get_gbox_p(geom2, &box2) )
2011	{
2012	if ( ( box2.xmin < box1.xmin ) \|\| ( box2.xmax > box1.xmax ) \|\|
2013	( box2.ymin < box1.ymin ) \|\| ( box2.ymax > box1.ymax ) )
2014	{
2015	PG_RETURN_BOOL(FALSE);
2016	}
2017	}
2018
2019	/*
2020	** short-circuit 2: if geom2 is a point and geom1 is a polygon
2021	** call the point-in-polygon function.
2022	*/
2023	type1 = gserialized_get_type(geom1);
2024	type2 = gserialized_get_type(geom2);
2025	if ((type1 == POLYGONTYPE \|\| type1 == MULTIPOLYGONTYPE) && type2 == POINTTYPE)
2026	{
2027	POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2028	lwgeom = lwgeom_from_gserialized(geom1);
2029	point = lwgeom_as_lwpoint(lwgeom_from_gserialized(geom2));
2030
2031	POSTGIS_DEBUGF(3, "Precall point_in_multipolygon_rtree %p, %p", lwgeom, point);
2032
2033	poly_cache = GetRtreeCache(fcinfo, lwgeom, geom1);
2034
2035	if ( poly_cache->ringIndices )
2036	{
2037	result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
2038	}
2039	else if ( type1 == POLYGONTYPE )
2040	{
2041	result = point_in_polygon((LWPOLY*)lwgeom, point);
2042	}
2043	else if ( type1 == MULTIPOLYGONTYPE )
2044	{
2045	result = point_in_multipolygon((LWMPOLY*)lwgeom, point);
2046	}
2047	else
2048	{
2049	/* Gulp! Should not be here... */
2050	elog(ERROR,"Type isn't poly or multipoly!");
2051	PG_RETURN_NULL();
2052	}
2053	lwgeom_free(lwgeom);
2054	lwpoint_free(point);
2055	PG_FREE_IF_COPY(geom1, 0);
2056	PG_FREE_IF_COPY(geom2, 1);
2057	if ( result == 1 ) /* completely inside */
2058	{
2059	PG_RETURN_BOOL(TRUE);
2060	}
2061	else
2062	{
2063	PG_RETURN_BOOL(FALSE);
2064	}
2065	}
2066	else
2067	{
2068	POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", type1, type2);
2069	}
2070
2071	initGEOS(lwnotice, lwgeom_geos_error);
2072
2073	#ifdef PREPARED_GEOM
2074	prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
2075
2076	if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
2077	{
2078	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2079	if ( 0 == g1 ) /* exception thrown at construction */
2080	{
2081	lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2082	PG_RETURN_NULL();
2083	}
2084	POSTGIS_DEBUG(4, "containsPrepared: cache is live, running preparedcontains");
2085	result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
2086	GEOSGeom_destroy(g1);
2087	}
2088	else
2089	#endif
2090	{
2091	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2092	if ( 0 == g1 ) /* exception thrown at construction */
2093	{
2094	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2095	PG_RETURN_NULL();
2096	}
2097	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2098	if ( 0 == g2 ) /* exception thrown at construction */
2099	{
2100	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2101	GEOSGeom_destroy(g1);
2102	PG_RETURN_NULL();
2103	}
2104	POSTGIS_DEBUG(4, "containsPrepared: cache is not ready, running standard contains");
2105	result = GEOSContains( g1, g2);
2106	GEOSGeom_destroy(g1);
2107	GEOSGeom_destroy(g2);
2108	}
2109
2110	if (result == 2)
2111	{
2112	lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
2113	PG_RETURN_NULL(); /* never get here */
2114	}
2115
2116	PG_FREE_IF_COPY(geom1, 0);
2117	PG_FREE_IF_COPY(geom2, 1);
2118
2119	PG_RETURN_BOOL(result);
2120
2121	}
2122
2123	PG_FUNCTION_INFO_V1(containsproperly);
2124	Datum containsproperly(PG_FUNCTION_ARGS)
2125	{
2126	GSERIALIZED * geom1;
2127	GSERIALIZED * geom2;
2128	bool result;
2129	GBOX box1, box2;
2130	#ifdef PREPARED_GEOM
2131	PrepGeomCache * prep_cache;
2132	#endif
2133
2134	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2135	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2136
2137	errorIfGeometryCollection(geom1,geom2);
2138	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2139
2140	/* A.ContainsProperly(Empty) == FALSE */
2141	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2142	PG_RETURN_BOOL(false);
2143
2144	/*
2145	* short-circuit: if geom2 bounding box is not completely inside
2146	* geom1 bounding box we can prematurely return FALSE.
2147	* Do the test IFF BOUNDING BOX AVAILABLE.
2148	*/
2149	if ( gserialized_get_gbox_p(geom1, &box1) &&
2150	gserialized_get_gbox_p(geom2, &box2) )
2151	{
2152	if (( box2.xmin < box1.xmin ) \|\| ( box2.xmax > box1.xmax ) \|\|
2153	( box2.ymin < box1.ymin ) \|\| ( box2.ymax > box1.ymax ))
2154	PG_RETURN_BOOL(FALSE);
2155	}
2156
2157	initGEOS(lwnotice, lwgeom_geos_error);
2158
2159	#ifdef PREPARED_GEOM
2160	prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
2161
2162	if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
2163	{
2164	GEOSGeometry g = (GEOSGeometry )POSTGIS2GEOS(geom2);
2165	if ( 0 == g ) /* exception thrown at construction */
2166	{
2167	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2168	PG_RETURN_NULL();
2169	}
2170	result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
2171	GEOSGeom_destroy(g);
2172	}
2173	else
2174	#endif
2175	{
2176	GEOSGeometry *g2;
2177	GEOSGeometry *g1;
2178
2179	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2180	if ( 0 == g1 ) /* exception thrown at construction */
2181	{
2182	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2183	PG_RETURN_NULL();
2184	}
2185	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2186	if ( 0 == g2 ) /* exception thrown at construction */
2187	{
2188	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2189	GEOSGeom_destroy(g1);
2190	PG_RETURN_NULL();
2191	}
2192	result = GEOSRelatePattern( g1, g2, "T*FFFF*" );
2193	GEOSGeom_destroy(g1);
2194	GEOSGeom_destroy(g2);
2195	}
2196
2197	if (result == 2)
2198	{
2199	lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
2200	PG_RETURN_NULL(); /* never get here */
2201	}
2202
2203	PG_FREE_IF_COPY(geom1, 0);
2204	PG_FREE_IF_COPY(geom2, 1);
2205
2206	PG_RETURN_BOOL(result);
2207	}
2208
2209	/*
2210	* Described at
2211	* http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2212	*/
2213	PG_FUNCTION_INFO_V1(covers);
2214	Datum covers(PG_FUNCTION_ARGS)
2215	{
2216	GSERIALIZED *geom1;
2217	GSERIALIZED *geom2;
2218	bool result;
2219	GBOX box1, box2;
2220	int type1, type2;
2221	LWGEOM *lwgeom;
2222	LWPOINT *point;
2223	RTREE_POLY_CACHE *poly_cache;
2224	#ifdef PREPARED_GEOM
2225	PrepGeomCache *prep_cache;
2226	#endif
2227
2228	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2229	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2230
2231	/* A.Covers(Empty) == FALSE */
2232	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2233	PG_RETURN_BOOL(false);
2234
2235	errorIfGeometryCollection(geom1,geom2);
2236	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2237
2238	/*
2239	* short-circuit 1: if geom2 bounding box is not completely inside
2240	* geom1 bounding box we can prematurely return FALSE.
2241	* Do the test IFF BOUNDING BOX AVAILABLE.
2242	*/
2243	if ( gserialized_get_gbox_p(geom1, &box1) &&
2244	gserialized_get_gbox_p(geom2, &box2) )
2245	{
2246	if (( box2.xmin < box1.xmin ) \|\| ( box2.xmax > box1.xmax ) \|\|
2247	( box2.ymin < box1.ymin ) \|\| ( box2.ymax > box1.ymax ))
2248	{
2249	PG_RETURN_BOOL(FALSE);
2250	}
2251	}
2252	/*
2253	* short-circuit 2: if geom2 is a point and geom1 is a polygon
2254	* call the point-in-polygon function.
2255	*/
2256	type1 = gserialized_get_type(geom1);
2257	type2 = gserialized_get_type(geom2);
2258	if ((type1 == POLYGONTYPE \|\| type1 == MULTIPOLYGONTYPE) && type2 == POINTTYPE)
2259	{
2260	POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2261
2262	lwgeom = lwgeom_from_gserialized(geom1);
2263	point = lwgeom_as_lwpoint(lwgeom_from_gserialized(geom2));
2264
2265	POSTGIS_DEBUGF(3, "Precall point_in_multipolygon_rtree %p, %p", lwgeom, point);
2266
2267	poly_cache = GetRtreeCache(fcinfo, lwgeom, geom1);
2268
2269	if ( poly_cache->ringIndices )
2270	{
2271	result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
2272	}
2273	else if ( type1 == POLYGONTYPE )
2274	{
2275	result = point_in_polygon((LWPOLY*)lwgeom, point);
2276	}
2277	else if ( type1 == MULTIPOLYGONTYPE )
2278	{
2279	result = point_in_multipolygon((LWMPOLY*)lwgeom, point);
2280	}
2281	else
2282	{
2283	/* Gulp! Should not be here... */
2284	elog(ERROR,"Type isn't poly or multipoly!");
2285	PG_RETURN_NULL();
2286	}
2287
2288	lwgeom_free(lwgeom);
2289	lwpoint_free(point);
2290	PG_FREE_IF_COPY(geom1, 0);
2291	PG_FREE_IF_COPY(geom2, 1);
2292	if ( result != -1 ) /* not outside */
2293	{
2294	PG_RETURN_BOOL(TRUE);
2295	}
2296	else
2297	{
2298	PG_RETURN_BOOL(FALSE);
2299	}
2300	}
2301	else
2302	{
2303	POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", type1, type2);
2304	}
2305
2306	initGEOS(lwnotice, lwgeom_geos_error);
2307
2308	#ifdef PREPARED_GEOM
2309	prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
2310
2311	if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
2312	{
2313	GEOSGeometry g1 = (GEOSGeometry )POSTGIS2GEOS(geom2);
2314	if ( 0 == g1 ) /* exception thrown at construction */
2315	{
2316	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2317	PG_RETURN_NULL();
2318	}
2319	result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
2320	GEOSGeom_destroy(g1);
2321	}
2322	else
2323	#endif
2324	{
2325	GEOSGeometry *g1;
2326	GEOSGeometry *g2;
2327
2328	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2329	if ( 0 == g1 ) /* exception thrown at construction */
2330	{
2331	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2332	PG_RETURN_NULL();
2333	}
2334	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2335	if ( 0 == g2 ) /* exception thrown at construction */
2336	{
2337	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2338	GEOSGeom_destroy(g1);
2339	PG_RETURN_NULL();
2340	}
2341	result = GEOSRelatePattern( g1, g2, "*****FF" );
2342	GEOSGeom_destroy(g1);
2343	GEOSGeom_destroy(g2);
2344	}
2345
2346	if (result == 2)
2347	{
2348	lwerror("GEOSCovers: %s", lwgeom_geos_errmsg);
2349	PG_RETURN_NULL(); /* never get here */
2350	}
2351
2352	PG_FREE_IF_COPY(geom1, 0);
2353	PG_FREE_IF_COPY(geom2, 1);
2354
2355	PG_RETURN_BOOL(result);
2356
2357	}
2358
2359
2360	/**
2361	* ST_Within(A, B) => ST_Contains(B, A) so we just delegate this calculation to the
2362	* Contains implementation.
2363	PG_FUNCTION_INFO_V1(within);
2364	Datum within(PG_FUNCTION_ARGS)
2365	*/
2366
2367	/*
2368	* Described at:
2369	* http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2370	*/
2371	PG_FUNCTION_INFO_V1(coveredby);
2372	Datum coveredby(PG_FUNCTION_ARGS)
2373	{
2374	GSERIALIZED *geom1;
2375	GSERIALIZED *geom2;
2376	GEOSGeometry g1, g2;
2377	bool result;
2378	GBOX box1, box2;
2379	LWGEOM *lwgeom;
2380	LWPOINT *point;
2381	int type1, type2;
2382	RTREE_POLY_CACHE *poly_cache;
2383	char patt = "FF**";
2384
2385	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2386	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2387
2388	errorIfGeometryCollection(geom1,geom2);
2389	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2390
2391	/* A.CoveredBy(Empty) == FALSE */
2392	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2393	PG_RETURN_BOOL(false);
2394
2395	/*
2396	* short-circuit 1: if geom1 bounding box is not completely inside
2397	* geom2 bounding box we can prematurely return FALSE.
2398	* Do the test IFF BOUNDING BOX AVAILABLE.
2399	*/
2400	if ( gserialized_get_gbox_p(geom1, &box1) &&
2401	gserialized_get_gbox_p(geom2, &box2) )
2402	{
2403	if ( ( box1.xmin < box2.xmin ) \|\| ( box1.xmax > box2.xmax ) \|\|
2404	( box1.ymin < box2.ymin ) \|\| ( box1.ymax > box2.ymax ) )
2405	{
2406	PG_RETURN_BOOL(FALSE);
2407	}
2408
2409	POSTGIS_DEBUG(3, "bounding box short-circuit missed.");
2410	}
2411	/*
2412	* short-circuit 2: if geom1 is a point and geom2 is a polygon
2413	* call the point-in-polygon function.
2414	*/
2415	type1 = gserialized_get_type(geom1);
2416	type2 = gserialized_get_type(geom2);
2417	if ((type2 == POLYGONTYPE \|\| type2 == MULTIPOLYGONTYPE) && type1 == POINTTYPE)
2418	{
2419	POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2420
2421	point = lwgeom_as_lwpoint(lwgeom_from_gserialized(geom1));
2422	lwgeom = lwgeom_from_gserialized(geom2);
2423
2424	poly_cache = GetRtreeCache(fcinfo, lwgeom, geom2);
2425
2426	if ( poly_cache->ringIndices )
2427	{
2428	result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
2429	}
2430	else if ( type2 == POLYGONTYPE )
2431	{
2432	result = point_in_polygon((LWPOLY*)lwgeom, point);
2433	}
2434	else if ( type2 == MULTIPOLYGONTYPE )
2435	{
2436	result = point_in_multipolygon((LWMPOLY*)lwgeom, point);
2437	}
2438	else
2439	{
2440	/* Gulp! Should not be here... */
2441	elog(ERROR,"Type isn't poly or multipoly!");
2442	PG_RETURN_NULL();
2443	}
2444
2445	lwgeom_free(lwgeom);
2446	lwpoint_free(point);
2447	PG_FREE_IF_COPY(geom1, 0);
2448	PG_FREE_IF_COPY(geom2, 1);
2449	if ( result != -1 ) /* not outside */
2450	{
2451	PG_RETURN_BOOL(TRUE);
2452	}
2453	else
2454	{
2455	PG_RETURN_BOOL(FALSE);
2456	}
2457	}
2458
2459	initGEOS(lwnotice, lwgeom_geos_error);
2460
2461	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2462
2463	if ( 0 == g1 ) /* exception thrown at construction */
2464	{
2465	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2466	PG_RETURN_NULL();
2467	}
2468
2469	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2470
2471	if ( 0 == g2 ) /* exception thrown at construction */
2472	{
2473	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2474	GEOSGeom_destroy(g1);
2475	PG_RETURN_NULL();
2476	}
2477
2478	result = GEOSRelatePattern(g1,g2,patt);
2479
2480	GEOSGeom_destroy(g1);
2481	GEOSGeom_destroy(g2);
2482
2483	if (result == 2)
2484	{
2485	lwerror("GEOSCoveredBy: %s", lwgeom_geos_errmsg);
2486	PG_RETURN_NULL(); /* never get here */
2487	}
2488
2489	PG_FREE_IF_COPY(geom1, 0);
2490	PG_FREE_IF_COPY(geom2, 1);
2491
2492	PG_RETURN_BOOL(result);
2493	}
2494
2495
2496
2497	PG_FUNCTION_INFO_V1(crosses);
2498	Datum crosses(PG_FUNCTION_ARGS)
2499	{
2500	GSERIALIZED *geom1;
2501	GSERIALIZED *geom2;
2502	GEOSGeometry g1, g2;
2503	bool result;
2504	GBOX box1, box2;
2505
2506	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2507	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2508
2509	errorIfGeometryCollection(geom1,geom2);
2510	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2511
2512	/* A.Crosses(Empty) == FALSE */
2513	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2514	PG_RETURN_BOOL(false);
2515
2516	/*
2517	* short-circuit 1: if geom2 bounding box does not overlap
2518	* geom1 bounding box we can prematurely return FALSE.
2519	* Do the test IFF BOUNDING BOX AVAILABLE.
2520	*/
2521	if ( gserialized_get_gbox_p(geom1, &box1) &&
2522	gserialized_get_gbox_p(geom2, &box2) )
2523	{
2524	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2525	{
2526	PG_RETURN_BOOL(FALSE);
2527	}
2528	}
2529
2530	initGEOS(lwnotice, lwgeom_geos_error);
2531
2532	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2533	if ( 0 == g1 ) /* exception thrown at construction */
2534	{
2535	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2536	PG_RETURN_NULL();
2537	}
2538
2539	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2540	if ( 0 == g2 ) /* exception thrown at construction */
2541	{
2542	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2543	GEOSGeom_destroy(g1);
2544	PG_RETURN_NULL();
2545	}
2546
2547	result = GEOSCrosses(g1,g2);
2548
2549	GEOSGeom_destroy(g1);
2550	GEOSGeom_destroy(g2);
2551
2552	if (result == 2)
2553	{
2554	lwerror("GEOSCrosses: %s", lwgeom_geos_errmsg);
2555	PG_RETURN_NULL(); /* never get here */
2556	}
2557
2558	PG_FREE_IF_COPY(geom1, 0);
2559	PG_FREE_IF_COPY(geom2, 1);
2560
2561	PG_RETURN_BOOL(result);
2562	}
2563
2564
2565	PG_FUNCTION_INFO_V1(intersects);
2566	Datum intersects(PG_FUNCTION_ARGS)
2567	{
2568	GSERIALIZED *geom1;
2569	GSERIALIZED *geom2;
2570	GSERIALIZED *serialized_poly;
2571	bool result;
2572	GBOX box1, box2;
2573	int type1, type2, polytype;
2574	LWPOINT *point;
2575	LWGEOM *lwgeom;
2576	RTREE_POLY_CACHE *poly_cache;
2577	#ifdef PREPARED_GEOM
2578	PrepGeomCache *prep_cache;
2579	#endif
2580
2581	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2582	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2583
2584	errorIfGeometryCollection(geom1,geom2);
2585	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2586
2587	/* A.Intersects(Empty) == FALSE */
2588	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2589	PG_RETURN_BOOL(false);
2590
2591	/*
2592	* short-circuit 1: if geom2 bounding box does not overlap
2593	* geom1 bounding box we can prematurely return FALSE.
2594	* Do the test IFF BOUNDING BOX AVAILABLE.
2595	*/
2596	if ( gserialized_get_gbox_p(geom1, &box1) &&
2597	gserialized_get_gbox_p(geom2, &box2) )
2598	{
2599	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2600	{
2601	PG_RETURN_BOOL(FALSE);
2602	}
2603	}
2604
2605	/*
2606	* short-circuit 2: if the geoms are a point and a polygon,
2607	* call the point_outside_polygon function.
2608	*/
2609	type1 = gserialized_get_type(geom1);
2610	type2 = gserialized_get_type(geom2);
2611	if ( (type1 == POINTTYPE && (type2 == POLYGONTYPE \|\| type2 == MULTIPOLYGONTYPE)) \|\|
2612	(type2 == POINTTYPE && (type1 == POLYGONTYPE \|\| type1 == MULTIPOLYGONTYPE)))
2613	{
2614	POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2615
2616	if ( type1 == POINTTYPE )
2617	{
2618	point = lwgeom_as_lwpoint(lwgeom_from_gserialized(geom1));
2619	lwgeom = lwgeom_from_gserialized(geom2);
2620	serialized_poly = geom2;
2621	polytype = type2;
2622	}
2623	else
2624	{
2625	point = lwgeom_as_lwpoint(lwgeom_from_gserialized(geom2));
2626	lwgeom = lwgeom_from_gserialized(geom1);
2627	serialized_poly = geom1;
2628	polytype = type1;
2629	}
2630
2631	poly_cache = GetRtreeCache(fcinfo, lwgeom, serialized_poly);
2632
2633	if ( poly_cache->ringIndices )
2634	{
2635	result = point_in_multipolygon_rtree(poly_cache->ringIndices, poly_cache->polyCount, poly_cache->ringCounts, point);
2636	}
2637	else if ( polytype == POLYGONTYPE )
2638	{
2639	result = point_in_polygon((LWPOLY*)lwgeom, point);
2640	}
2641	else if ( polytype == MULTIPOLYGONTYPE )
2642	{
2643	result = point_in_multipolygon((LWMPOLY*)lwgeom, point);
2644	}
2645	else
2646	{
2647	/* Gulp! Should not be here... */
2648	elog(ERROR,"Type isn't poly or multipoly!");
2649	PG_RETURN_NULL();
2650	}
2651
2652	lwgeom_free(lwgeom);
2653	lwpoint_free(point);
2654	PG_FREE_IF_COPY(geom1, 0);
2655	PG_FREE_IF_COPY(geom2, 1);
2656	if ( result != -1 ) /* not outside */
2657	{
2658	PG_RETURN_BOOL(TRUE);
2659	}
2660	else
2661	{
2662	PG_RETURN_BOOL(FALSE);
2663	}
2664	}
2665
2666	initGEOS(lwnotice, lwgeom_geos_error);
2667	#ifdef PREPARED_GEOM
2668	prep_cache = GetPrepGeomCache( fcinfo, geom1, geom2 );
2669
2670	if ( prep_cache && prep_cache->prepared_geom )
2671	{
2672	if ( prep_cache->argnum == 1 )
2673	{
2674	GEOSGeometry g = (GEOSGeometry )POSTGIS2GEOS(geom2);
2675	if ( 0 == g ) /* exception thrown at construction */
2676	{
2677	lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2678	PG_RETURN_NULL();
2679	}
2680	result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2681	GEOSGeom_destroy(g);
2682	}
2683	else
2684	{
2685	GEOSGeometry g = (GEOSGeometry )POSTGIS2GEOS(geom1);
2686	if ( 0 == g ) /* exception thrown at construction */
2687	{
2688	lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2689	PG_RETURN_NULL();
2690	}
2691	result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2692	GEOSGeom_destroy(g);
2693	}
2694	}
2695	else
2696	#endif
2697	{
2698	GEOSGeometry *g1;
2699	GEOSGeometry *g2;
2700	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2701	if ( 0 == g1 ) /* exception thrown at construction */
2702	{
2703	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2704	PG_RETURN_NULL();
2705	}
2706	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2707	if ( 0 == g2 ) /* exception thrown at construction */
2708	{
2709	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2710	GEOSGeom_destroy(g1);
2711	PG_RETURN_NULL();
2712	}
2713	result = GEOSIntersects( g1, g2);
2714	GEOSGeom_destroy(g1);
2715	GEOSGeom_destroy(g2);
2716	}
2717
2718	if (result == 2)
2719	{
2720	lwerror("GEOSIntersects: %s", lwgeom_geos_errmsg);
2721	PG_RETURN_NULL(); /* never get here */
2722	}
2723
2724	PG_FREE_IF_COPY(geom1, 0);
2725	PG_FREE_IF_COPY(geom2, 1);
2726
2727	PG_RETURN_BOOL(result);
2728	}
2729
2730
2731	PG_FUNCTION_INFO_V1(touches);
2732	Datum touches(PG_FUNCTION_ARGS)
2733	{
2734	GSERIALIZED *geom1;
2735	GSERIALIZED *geom2;
2736	GEOSGeometry g1, g2;
2737	bool result;
2738	GBOX box1, box2;
2739
2740	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2741	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2742
2743	errorIfGeometryCollection(geom1,geom2);
2744	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2745
2746	/* A.Touches(Empty) == FALSE */
2747	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2748	PG_RETURN_BOOL(false);
2749
2750	/*
2751	* short-circuit 1: if geom2 bounding box does not overlap
2752	* geom1 bounding box we can prematurely return FALSE.
2753	* Do the test IFF BOUNDING BOX AVAILABLE.
2754	*/
2755	if ( gserialized_get_gbox_p(geom1, &box1) &&
2756	gserialized_get_gbox_p(geom2, &box2) )
2757	{
2758	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2759	{
2760	PG_RETURN_BOOL(FALSE);
2761	}
2762	}
2763
2764	initGEOS(lwnotice, lwgeom_geos_error);
2765
2766	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1 );
2767	if ( 0 == g1 ) /* exception thrown at construction */
2768	{
2769	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2770	PG_RETURN_NULL();
2771	}
2772
2773	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2 );
2774	if ( 0 == g2 ) /* exception thrown at construction */
2775	{
2776	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2777	GEOSGeom_destroy(g1);
2778	PG_RETURN_NULL();
2779	}
2780
2781	result = GEOSTouches(g1,g2);
2782
2783	GEOSGeom_destroy(g1);
2784	GEOSGeom_destroy(g2);
2785
2786	if (result == 2)
2787	{
2788	lwerror("GEOSTouches: %s", lwgeom_geos_errmsg);
2789	PG_RETURN_NULL(); /* never get here */
2790	}
2791
2792	PG_FREE_IF_COPY(geom1, 0);
2793	PG_FREE_IF_COPY(geom2, 1);
2794
2795	PG_RETURN_BOOL(result);
2796	}
2797
2798
2799	PG_FUNCTION_INFO_V1(disjoint);
2800	Datum disjoint(PG_FUNCTION_ARGS)
2801	{
2802	GSERIALIZED *geom1;
2803	GSERIALIZED *geom2;
2804	GEOSGeometry g1, g2;
2805	bool result;
2806	GBOX box1, box2;
2807
2808	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2809	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2810
2811	errorIfGeometryCollection(geom1,geom2);
2812	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2813
2814	/* A.Disjoint(Empty) == TRUE */
2815	if ( gserialized_is_empty(geom1) \|\| gserialized_is_empty(geom2) )
2816	PG_RETURN_BOOL(true);
2817
2818	/*
2819	* short-circuit 1: if geom2 bounding box does not overlap
2820	* geom1 bounding box we can prematurely return TRUE.
2821	* Do the test IFF BOUNDING BOX AVAILABLE.
2822	*/
2823	if ( gserialized_get_gbox_p(geom1, &box1) &&
2824	gserialized_get_gbox_p(geom2, &box2) )
2825	{
2826	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2827	{
2828	PG_RETURN_BOOL(TRUE);
2829	}
2830	}
2831
2832	initGEOS(lwnotice, lwgeom_geos_error);
2833
2834	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2835	if ( 0 == g1 ) /* exception thrown at construction */
2836	{
2837	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2838	PG_RETURN_NULL();
2839	}
2840
2841	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2842	if ( 0 == g2 ) /* exception thrown at construction */
2843	{
2844	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2845	GEOSGeom_destroy(g1);
2846	PG_RETURN_NULL();
2847	}
2848
2849	result = GEOSDisjoint(g1,g2);
2850
2851	GEOSGeom_destroy(g1);
2852	GEOSGeom_destroy(g2);
2853
2854	if (result == 2)
2855	{
2856	lwerror("GEOSDisjoint: %s", lwgeom_geos_errmsg);
2857	PG_RETURN_NULL(); /* never get here */
2858	}
2859
2860	PG_FREE_IF_COPY(geom1, 0);
2861	PG_FREE_IF_COPY(geom2, 1);
2862
2863	PG_RETURN_BOOL(result);
2864	}
2865
2866
2867	PG_FUNCTION_INFO_V1(relate_pattern);
2868	Datum relate_pattern(PG_FUNCTION_ARGS)
2869	{
2870	GSERIALIZED *geom1;
2871	GSERIALIZED *geom2;
2872	char *patt;
2873	bool result;
2874	GEOSGeometry g1, g2;
2875	int i;
2876
2877	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2878	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2879
2880
2881	/* TODO handle empty */
2882
2883	errorIfGeometryCollection(geom1,geom2);
2884	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2885
2886	initGEOS(lwnotice, lwgeom_geos_error);
2887
2888	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
2889	if ( 0 == g1 ) /* exception thrown at construction */
2890	{
2891	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2892	PG_RETURN_NULL();
2893	}
2894	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
2895	if ( 0 == g2 ) /* exception thrown at construction */
2896	{
2897	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2898	GEOSGeom_destroy(g1);
2899	PG_RETURN_NULL();
2900	}
2901
2902	patt = DatumGetCString(DirectFunctionCall1(textout,
2903	PointerGetDatum(PG_GETARG_DATUM(2))));
2904
2905	/*
2906	** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2907	*/
2908	for ( i = 0; i < strlen(patt); i++ )
2909	{
2910	if ( patt[i] == 't' ) patt[i] = 'T';
2911	if ( patt[i] == 'f' ) patt[i] = 'F';
2912	}
2913
2914	result = GEOSRelatePattern(g1,g2,patt);
2915	GEOSGeom_destroy(g1);
2916	GEOSGeom_destroy(g2);
2917	pfree(patt);
2918
2919	if (result == 2)
2920	{
2921	lwerror("GEOSRelatePattern: %s", lwgeom_geos_errmsg);
2922	PG_RETURN_NULL(); /* never get here */
2923	}
2924
2925	PG_FREE_IF_COPY(geom1, 0);
2926	PG_FREE_IF_COPY(geom2, 1);
2927
2928	PG_RETURN_BOOL(result);
2929	}
2930
2931
2932
2933	PG_FUNCTION_INFO_V1(relate_full);
2934	Datum relate_full(PG_FUNCTION_ARGS)
2935	{
2936	GSERIALIZED *geom1;
2937	GSERIALIZED *geom2;
2938	GEOSGeometry g1, g2;
2939	char *relate_str;
2940	text *result;
2941	#if POSTGIS_GEOS_VERSION >= 33
2942	int bnr = GEOSRELATE_BNR_OGC;
2943	#endif
2944
2945	POSTGIS_DEBUG(2, "in relate_full()");
2946
2947	/* TODO handle empty */
2948
2949	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2950	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2951
2952	if ( PG_NARGS() > 2 ) {
2953	#if POSTGIS_GEOS_VERSION >= 33
2954	bnr = PG_GETARG_INT32(2);
2955	#else
2956	lwerror("The GEOS version this PostGIS binary "
2957	"was compiled against (%d) doesn't support "
2958	"specifying a boundary node rule with ST_Relate"
2959	" (3.3.0+ required)",
2960	POSTGIS_GEOS_VERSION);
2961	PG_RETURN_NULL();
2962	#endif
2963	}
2964
2965	errorIfGeometryCollection(geom1,geom2);
2966	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
2967
2968	initGEOS(lwnotice, lwgeom_geos_error);
2969
2970	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1 );
2971	if ( 0 == g1 ) /* exception thrown at construction */
2972	{
2973	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2974	PG_RETURN_NULL();
2975	}
2976	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2 );
2977	if ( 0 == g2 ) /* exception thrown at construction */
2978	{
2979	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
2980	GEOSGeom_destroy(g1);
2981	PG_RETURN_NULL();
2982	}
2983
2984	POSTGIS_DEBUG(3, "constructed geometries ");
2985
2986	if ((g1==NULL) \|\| (g2 == NULL))
2987	elog(NOTICE,"g1 or g2 are null");
2988
2989	POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2990	POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2991
2992	#if POSTGIS_GEOS_VERSION >= 33
2993	relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2994	#else
2995	relate_str = GEOSRelate(g1, g2);
2996	#endif
2997
2998	GEOSGeom_destroy(g1);
2999	GEOSGeom_destroy(g2);
3000
3001	if (relate_str == NULL)
3002	{
3003	lwerror("GEOSRelate: %s", lwgeom_geos_errmsg);
3004	PG_RETURN_NULL(); /* never get here */
3005	}
3006
3007	result = cstring2text(relate_str);
3008	GEOSFree(relate_str);
3009
3010	PG_FREE_IF_COPY(geom1, 0);
3011	PG_FREE_IF_COPY(geom2, 1);
3012
3013	PG_RETURN_TEXT_P(result);
3014	}
3015
3016
3017	PG_FUNCTION_INFO_V1(ST_Equals);
3018	Datum ST_Equals(PG_FUNCTION_ARGS)
3019	{
3020	GSERIALIZED *geom1;
3021	GSERIALIZED *geom2;
3022	GEOSGeometry g1, g2;
3023	bool result;
3024	GBOX box1, box2;
3025
3026	geom1 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3027	geom2 = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3028
3029	errorIfGeometryCollection(geom1,geom2);
3030	error_if_srid_mismatch(gserialized_get_srid(geom1), gserialized_get_srid(geom2));
3031
3032	/* Empty == Empty */
3033	if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
3034	PG_RETURN_BOOL(TRUE);
3035
3036	/*
3037	* short-circuit: Loose test, if geom2 bounding box does not overlap
3038	* geom1 bounding box we can prematurely return FALSE.
3039	*
3040	* TODO: use gbox_same_2d instead (not available at time of writing)
3041	*/
3042	if ( gserialized_get_gbox_p(geom1, &box1) &&
3043	gserialized_get_gbox_p(geom2, &box2) )
3044	{
3045	if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
3046	{
3047	PG_RETURN_BOOL(FALSE);
3048	}
3049	}
3050
3051	initGEOS(lwnotice, lwgeom_geos_error);
3052
3053	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
3054
3055	if ( 0 == g1 ) /* exception thrown at construction */
3056	{
3057	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3058	PG_RETURN_NULL();
3059	}
3060
3061	g2 = (GEOSGeometry *)POSTGIS2GEOS(geom2);
3062
3063	if ( 0 == g2 ) /* exception thrown at construction */
3064	{
3065	lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3066	GEOSGeom_destroy(g1);
3067	PG_RETURN_NULL();
3068	}
3069
3070	result = GEOSEquals(g1,g2);
3071
3072	GEOSGeom_destroy(g1);
3073	GEOSGeom_destroy(g2);
3074
3075	if (result == 2)
3076	{
3077	lwerror("GEOSEquals: %s", lwgeom_geos_errmsg);
3078	PG_RETURN_NULL(); /never get here /
3079	}
3080
3081	PG_FREE_IF_COPY(geom1, 0);
3082	PG_FREE_IF_COPY(geom2, 1);
3083
3084	PG_RETURN_BOOL(result);
3085	}
3086
3087	PG_FUNCTION_INFO_V1(issimple);
3088	Datum issimple(PG_FUNCTION_ARGS)
3089	{
3090	GSERIALIZED *geom;
3091	GEOSGeometry *g1;
3092	int result;
3093
3094	POSTGIS_DEBUG(2, "issimple called");
3095
3096	geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3097
3098	if ( gserialized_is_empty(geom) )
3099	PG_RETURN_BOOL(TRUE);
3100
3101	initGEOS(lwnotice, lwgeom_geos_error);
3102
3103	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom);
3104	if ( 0 == g1 ) /* exception thrown at construction */
3105	{
3106	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3107	PG_RETURN_NULL();
3108	}
3109	result = GEOSisSimple(g1);
3110	GEOSGeom_destroy(g1);
3111
3112	if (result == 2)
3113	{
3114	lwerror("GEOSisSimple: %s", lwgeom_geos_errmsg);
3115	PG_RETURN_NULL(); /never get here /
3116	}
3117
3118	PG_FREE_IF_COPY(geom, 0);
3119
3120	PG_RETURN_BOOL(result);
3121	}
3122
3123	PG_FUNCTION_INFO_V1(isring);
3124	Datum isring(PG_FUNCTION_ARGS)
3125	{
3126	GSERIALIZED *geom;
3127	GEOSGeometry *g1;
3128	int result;
3129
3130	geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3131
3132	if (gserialized_get_type(geom) != LINETYPE)
3133	{
3134	elog(ERROR,"isring() should only be called on a LINE");
3135	}
3136
3137	/* Empty things can't close */
3138	if ( gserialized_is_empty(geom) )
3139	PG_RETURN_BOOL(FALSE);
3140
3141	initGEOS(lwnotice, lwgeom_geos_error);
3142
3143	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom );
3144	if ( 0 == g1 ) /* exception thrown at construction */
3145	{
3146	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3147	PG_RETURN_NULL();
3148	}
3149	result = GEOSisRing(g1);
3150	GEOSGeom_destroy(g1);
3151
3152	if (result == 2)
3153	{
3154	lwerror("GEOSisRing: %s", lwgeom_geos_errmsg);
3155	PG_RETURN_NULL();
3156	}
3157
3158	PG_FREE_IF_COPY(geom, 0);
3159	PG_RETURN_BOOL(result);
3160	}
3161
3162
3163
3164
3165
3166	GSERIALIZED *
3167	GEOS2POSTGIS(GEOSGeom geom, char want3d)
3168	{
3169	LWGEOM *lwgeom;
3170	GSERIALIZED *result;
3171
3172	lwgeom = GEOS2LWGEOM(geom, want3d);
3173	if ( ! lwgeom )
3174	{
3175	lwerror("GEOS2POSTGIS: GEOS2LWGEOM returned NULL");
3176	return NULL;
3177	}
3178
3179	POSTGIS_DEBUGF(4, "GEOS2POSTGIS: GEOS2LWGEOM returned a %s", lwgeom_summary(lwgeom, 0));
3180
3181	if ( lwgeom_needs_bbox(lwgeom) == LW_TRUE )
3182	{
3183	lwgeom_add_bbox(lwgeom);
3184	}
3185
3186	result = geometry_serialize(lwgeom);
3187	lwgeom_free(lwgeom);
3188
3189	return result;
3190	}
3191
3192	/-----=POSTGIS2GEOS= /
3193
3194
3195	GEOSGeometry *
3196	POSTGIS2GEOS(GSERIALIZED *pglwgeom)
3197	{
3198	GEOSGeometry *ret;
3199	LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
3200	if ( ! lwgeom )
3201	{
3202	lwerror("POSTGIS2GEOS: unable to deserialize input");
3203	return NULL;
3204	}
3205	ret = LWGEOM2GEOS(lwgeom);
3206	lwgeom_free(lwgeom);
3207	if ( ! ret )
3208	{
3209	/* lwerror("POSTGIS2GEOS conversion failed"); */
3210	return NULL;
3211	}
3212	return ret;
3213	}
3214
3215
3216	PG_FUNCTION_INFO_V1(GEOSnoop);
3217	Datum GEOSnoop(PG_FUNCTION_ARGS)
3218	{
3219	GSERIALIZED *geom;
3220	GEOSGeometry *geosgeom;
3221	GSERIALIZED *lwgeom_result;
3222	#if POSTGIS_DEBUG_LEVEL > 0
3223	int result;
3224	LWGEOM_UNPARSER_RESULT lwg_unparser_result;
3225	#endif
3226
3227	initGEOS(lwnotice, lwgeom_geos_error);
3228
3229	geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3230
3231
3232	geosgeom = (GEOSGeometry *)POSTGIS2GEOS(geom);
3233	if ( ! geosgeom ) PG_RETURN_NULL();
3234
3235	lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
3236	GEOSGeom_destroy(geosgeom);
3237
3238
3239	PG_FREE_IF_COPY(geom, 0);
3240
3241	PG_RETURN_POINTER(lwgeom_result);
3242	}
3243
3244	PG_FUNCTION_INFO_V1(polygonize_garray);
3245	Datum polygonize_garray(PG_FUNCTION_ARGS)
3246	{
3247	Datum datum;
3248	ArrayType *array;
3249	int is3d = 0;
3250	uint32 nelems, i;
3251	GSERIALIZED *result;
3252	GEOSGeometry *geos_result;
3253	const GEOSGeometry **vgeoms;
3254	int srid=SRID_UNKNOWN;
3255	size_t offset;
3256	#if POSTGIS_DEBUG_LEVEL >= 3
3257	static int call=1;
3258	#endif
3259
3260	#if POSTGIS_DEBUG_LEVEL >= 3
3261	call++;
3262	#endif
3263
3264	datum = PG_GETARG_DATUM(0);
3265
3266	/* Null array, null geometry (should be empty?) */
3267	if ( (Pointer *)datum == NULL ) PG_RETURN_NULL();
3268
3269	array = DatumGetArrayTypeP(datum);
3270
3271	nelems = ArrayGetNItems(ARR_NDIM(array), ARR_DIMS(array));
3272
3273	POSTGIS_DEBUGF(3, "polygonize_garray: number of elements: %d", nelems);
3274
3275	if ( nelems == 0 ) PG_RETURN_NULL();
3276
3277	/* Ok, we really need geos now ;) */
3278	initGEOS(lwnotice, lwgeom_geos_error);
3279
3280	vgeoms = palloc(sizeof(GEOSGeometry )nelems);
3281	offset = 0;
3282	for (i=0; i<nelems; i++)
3283	{
3284	GEOSGeometry* g;
3285	GSERIALIZED geom = (GSERIALIZED )(ARR_DATA_PTR(array)+offset);
3286	offset += INTALIGN(VARSIZE(geom));
3287	if ( ! is3d ) is3d = gserialized_has_z(geom);
3288
3289	g = (GEOSGeometry *)POSTGIS2GEOS(geom);
3290	if ( 0 == g ) /* exception thrown at construction */
3291	{
3292	lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3293	PG_RETURN_NULL();
3294	}
3295	vgeoms[i] = g;
3296	if ( ! i )
3297	{
3298	srid = gserialized_get_srid(geom);
3299	}
3300	else
3301	{
3302	if ( srid != gserialized_get_srid(geom) )
3303	{
3304	elog(ERROR, "polygonize: operation on mixed SRID geometries");
3305	PG_RETURN_NULL();
3306	}
3307	}
3308	}
3309
3310	POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
3311
3312	geos_result = GEOSPolygonize(vgeoms, nelems);
3313
3314	POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
3315
3316	for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
3317	pfree(vgeoms);
3318
3319	if ( ! geos_result ) PG_RETURN_NULL();
3320
3321	GEOSSetSRID(geos_result, srid);
3322	result = GEOS2POSTGIS(geos_result, is3d);
3323	GEOSGeom_destroy(geos_result);
3324	if ( result == NULL )
3325	{
3326	elog(ERROR, "GEOS2POSTGIS returned an error");
3327	PG_RETURN_NULL(); /never get here /
3328	}
3329
3330	/compressType(result); /
3331
3332	PG_RETURN_POINTER(result);
3333
3334	}
3335
3336	PG_FUNCTION_INFO_V1(linemerge);
3337	Datum linemerge(PG_FUNCTION_ARGS)
3338	{
3339	GSERIALIZED *geom1;
3340	GEOSGeometry g1, g3;
3341	GSERIALIZED *result;
3342
3343	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3344
3345	initGEOS(lwnotice, lwgeom_geos_error);
3346
3347	g1 = (GEOSGeometry *)POSTGIS2GEOS(geom1);
3348
3349	if ( 0 == g1 ) /* exception thrown at construction */
3350	{
3351	lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg);
3352	PG_RETURN_NULL();
3353	}
3354
3355	g3 = GEOSLineMerge(g1);
3356
3357	if (g3 == NULL)
3358	{
3359	elog(ERROR,"GEOS LineMerge() threw an error!");
3360	GEOSGeom_destroy(g1);
3361	PG_RETURN_NULL(); /never get here /
3362	}
3363
3364
3365	POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
3366
3367	GEOSSetSRID(g3, gserialized_get_srid(geom1));
3368
3369	result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
3370
3371	if (result == NULL)
3372	{
3373	GEOSGeom_destroy(g1);
3374	GEOSGeom_destroy(g3);
3375	elog(ERROR,"GEOS LineMerge() threw an error (result postgis geometry formation)!");
3376	PG_RETURN_NULL(); /never get here /
3377	}
3378	GEOSGeom_destroy(g1);
3379	GEOSGeom_destroy(g3);
3380
3381
3382	/* compressType(result); */
3383
3384	PG_FREE_IF_COPY(geom1, 0);
3385
3386	PG_RETURN_POINTER(result);
3387	}
3388
3389	/*
3390	* Take a geometry and return an areal geometry
3391	* (Polygon or MultiPolygon).
3392	* Actually a wrapper around GEOSpolygonize,
3393	* transforming the resulting collection into
3394	* a valid polygon Geometry.
3395	*/
3396	PG_FUNCTION_INFO_V1(ST_BuildArea);
3397	Datum ST_BuildArea(PG_FUNCTION_ARGS)
3398	{
3399	GSERIALIZED *result;
3400	GSERIALIZED *geom;
3401	LWGEOM lwgeom_in, lwgeom_out;
3402
3403	geom = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3404	lwgeom_in = lwgeom_from_gserialized(geom);
3405
3406	lwgeom_out = lwgeom_buildarea(lwgeom_in);
3407	lwgeom_free(lwgeom_in) ;
3408
3409	if ( ! lwgeom_out ) {
3410	PG_FREE_IF_COPY(geom, 0);
3411	PG_RETURN_NULL();
3412	}
3413
3414	result = geometry_serialize(lwgeom_out) ;
3415	lwgeom_free(lwgeom_out) ;
3416
3417	PG_FREE_IF_COPY(geom, 0);
3418	PG_RETURN_POINTER(result);
3419	}
3420
3421	/*
3422	* ST_Snap
3423	*
3424	* Snap a geometry to another with a given tolerance
3425	*/
3426	Datum ST_Snap(PG_FUNCTION_ARGS);
3427	PG_FUNCTION_INFO_V1(ST_Snap);
3428	Datum ST_Snap(PG_FUNCTION_ARGS)
3429	{
3430	#if POSTGIS_GEOS_VERSION < 33
3431	lwerror("The GEOS version this PostGIS binary "
3432	"was compiled against (%d) doesn't support "
3433	"'ST_Snap' function (3.3.0+ required)",
3434	POSTGIS_GEOS_VERSION);
3435	PG_RETURN_NULL();
3436	#else /* POSTGIS_GEOS_VERSION >= 33 */
3437	GSERIALIZED geom1, geom2, *result;
3438	LWGEOM lwgeom1, lwgeom2, *lwresult;
3439	double tolerance;
3440
3441	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3442	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3443	tolerance = PG_GETARG_FLOAT8(2);
3444
3445	lwgeom1 = lwgeom_from_gserialized(geom1);
3446	lwgeom2 = lwgeom_from_gserialized(geom2);
3447
3448	lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3449	lwgeom_free(lwgeom1);
3450	lwgeom_free(lwgeom2);
3451	PG_FREE_IF_COPY(geom1, 0);
3452	PG_FREE_IF_COPY(geom2, 1);
3453
3454	result = geometry_serialize(lwresult);
3455	lwgeom_free(lwresult);
3456
3457	PG_RETURN_POINTER(result);
3458
3459	#endif /* POSTGIS_GEOS_VERSION >= 33 */
3460
3461	}
3462
3463	/*
3464	* ST_Split
3465	*
3466	* Split polygon by line, line by line, line by point.
3467	* Returns at most components as a collection.
3468	* First element of the collection is always the part which
3469	* remains after the cut, while the second element is the
3470	* part which has been cut out. We arbitrarely take the part
3471	* on the right of cut lines as the part which has been cut out.
3472	* For a line cut by a point the part which remains is the one
3473	* from start of the line to the cut point.
3474	*
3475	*
3476	* Author: Sandro Santilli <strk@keybit.net>
3477	*
3478	* Work done for Faunalia (http://www.faunalia.it) with fundings
3479	* from Regione Toscana - Sistema Informativo per il Governo
3480	* del Territorio e dell'Ambiente (RT-SIGTA).
3481	*
3482	* Thanks to the PostGIS community for sharing poly/line ideas [1]
3483	*
3484	* [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3485	*
3486	*/
3487	Datum ST_Split(PG_FUNCTION_ARGS);
3488	PG_FUNCTION_INFO_V1(ST_Split);
3489	Datum ST_Split(PG_FUNCTION_ARGS)
3490	{
3491	GSERIALIZED in, blade_in, *out;
3492	LWGEOM lwgeom_in, lwblade_in, *lwgeom_out;
3493
3494	in = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3495	lwgeom_in = lwgeom_from_gserialized(in);
3496
3497	blade_in = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3498	lwblade_in = lwgeom_from_gserialized(blade_in);
3499
3500	error_if_srid_mismatch(lwgeom_in->srid, lwblade_in->srid);
3501
3502	lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3503	lwgeom_free(lwgeom_in);
3504	lwgeom_free(lwblade_in);
3505
3506	if ( ! lwgeom_out )
3507	{
3508	PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3509	PG_FREE_IF_COPY(blade_in, 1);
3510	PG_RETURN_NULL();
3511	}
3512
3513	out = geometry_serialize(lwgeom_out);
3514	lwgeom_free(lwgeom_out);
3515	PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3516	PG_FREE_IF_COPY(blade_in, 1);
3517
3518	PG_RETURN_POINTER(out);
3519	}
3520
3521	/**********************************************************************
3522	*
3523	* ST_SharedPaths
3524	*
3525	* Return the set of paths shared between two linear geometries,
3526	* and their direction (same or opposite).
3527	*
3528	* Developed by Sandro Santilli (strk@keybit.net) for Faunalia
3529	* (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3530	* Informativo per la Gestione del Territorio e dell' Ambiente
3531	* [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3532	* trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3533	* 0494241492)"
3534	*
3535	**********************************************************************/
3536	Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3537	PG_FUNCTION_INFO_V1(ST_SharedPaths);
3538	Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3539	{
3540	#if POSTGIS_GEOS_VERSION < 33
3541	lwerror("The GEOS version this PostGIS binary "
3542	"was compiled against (%d) doesn't support "
3543	"'ST_SharedPaths' function (3.3.0+ required)",
3544	POSTGIS_GEOS_VERSION);
3545	PG_RETURN_NULL();
3546	#else /* POSTGIS_GEOS_VERSION >= 33 */
3547	GSERIALIZED geom1, geom2, *out;
3548	LWGEOM g1, g2, *lwgeom_out;
3549
3550	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3551	geom2 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
3552
3553	g1 = lwgeom_from_gserialized(geom1);
3554	g2 = lwgeom_from_gserialized(geom2);
3555
3556	lwgeom_out = lwgeom_sharedpaths(g1, g2);
3557	lwgeom_free(g1);
3558	lwgeom_free(g2);
3559
3560	if ( ! lwgeom_out )
3561	{
3562	PG_FREE_IF_COPY(geom1, 0);
3563	PG_FREE_IF_COPY(geom2, 1);
3564	PG_RETURN_NULL();
3565	}
3566
3567	out = geometry_serialize(lwgeom_out);
3568	lwgeom_free(lwgeom_out);
3569
3570	PG_FREE_IF_COPY(geom1, 0);
3571	PG_FREE_IF_COPY(geom2, 1);
3572	PG_RETURN_POINTER(out);
3573
3574	#endif /* POSTGIS_GEOS_VERSION >= 33 */
3575
3576	}
3577
3578
3579	/**********************************************************************
3580	*
3581	* ST_Node
3582	*
3583	* Fully node a set of lines using the least possible nodes while
3584	* preserving all of the input ones.
3585	*
3586	**********************************************************************/
3587	Datum ST_Node(PG_FUNCTION_ARGS);
3588	PG_FUNCTION_INFO_V1(ST_Node);
3589	Datum ST_Node(PG_FUNCTION_ARGS)
3590	{
3591	#if POSTGIS_GEOS_VERSION < 33
3592	lwerror("The GEOS version this PostGIS binary "
3593	"was compiled against (%d) doesn't support "
3594	"'ST_Node' function (3.3.0+ required)",
3595	POSTGIS_GEOS_VERSION);
3596	PG_RETURN_NULL();
3597	#else /* POSTGIS_GEOS_VERSION >= 33 */
3598	GSERIALIZED geom1, out;
3599	LWGEOM g1, lwgeom_out;
3600
3601	geom1 = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3602
3603	g1 = lwgeom_from_gserialized(geom1);
3604
3605	lwgeom_out = lwgeom_node(g1);
3606	lwgeom_free(g1);
3607
3608	if ( ! lwgeom_out )
3609	{
3610	PG_FREE_IF_COPY(geom1, 0);
3611	PG_RETURN_NULL();
3612	}
3613
3614	out = geometry_serialize(lwgeom_out);
3615	lwgeom_free(lwgeom_out);
3616
3617	PG_FREE_IF_COPY(geom1, 0);
3618	PG_RETURN_POINTER(out);
3619
3620	#endif /* POSTGIS_GEOS_VERSION >= 33 */
3621
3622	}

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format