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