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