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