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 strictly contained by bbox2, return input geometry */
1587  if (bbox2->xmin < bbox1.xmin && bbox2->xmax > bbox1.xmax &&
1588  bbox2->ymin < bbox1.ymin && bbox2->ymax > bbox1.ymax)
1589  {
1590  PG_RETURN_DATUM(PG_GETARG_DATUM(geom_idx));
1591  }
1592 
1593  /* If bbox1 outside of bbox2, return empty */
1594  if (!gbox_overlaps_2d(&bbox1, bbox2))
1595  {
1596  /* Get type and srid from datum */
1597  lwresult = lwgeom_construct_empty(type, srid, 0, 0);
1598  result = geometry_serialize(lwresult) ;
1599  lwgeom_free(lwresult) ;
1600  PG_RETURN_POINTER(result);
1601  }
1602 
1603  lwgeom1 = lwgeom_from_gserialized(PG_GETARG_GSERIALIZED_P(geom_idx));
1604  lwresult = lwgeom_clip_by_rect(lwgeom1, bbox2->xmin, bbox2->ymin,
1605  bbox2->xmax, bbox2->ymax);
1606 
1607  lwgeom_free(lwgeom1);
1608 
1609  if (!lwresult)
1610  PG_RETURN_NULL();
1611 
1612  result = geometry_serialize(lwresult) ;
1613  PG_RETURN_POINTER(result);
1614 }
1615 
1616 
1617 /*---------------------------------------------*/
1618 
1619 PG_FUNCTION_INFO_V1(isvalid);
1620 Datum isvalid(PG_FUNCTION_ARGS)
1621 {
1622  GSERIALIZED *geom1;
1623  LWGEOM *lwgeom;
1624  char result;
1625  GEOSGeom g1;
1626 
1627  geom1 = PG_GETARG_GSERIALIZED_P(0);
1628 
1629  /* Empty.IsValid() == TRUE */
1630  if ( gserialized_is_empty(geom1) )
1631  PG_RETURN_BOOL(true);
1632 
1633  initGEOS(lwpgnotice, lwgeom_geos_error);
1634 
1635  lwgeom = lwgeom_from_gserialized(geom1);
1636  if ( ! lwgeom )
1637  {
1638  lwpgerror("unable to deserialize input");
1639  }
1640  g1 = LWGEOM2GEOS(lwgeom, 0);
1641  lwgeom_free(lwgeom);
1642 
1643  if ( ! g1 )
1644  {
1645  /* should we drop the following
1646  * notice now that we have ST_isValidReason ?
1647  */
1648  lwpgnotice("%s", lwgeom_geos_errmsg);
1649  PG_RETURN_BOOL(false);
1650  }
1651 
1652  result = GEOSisValid(g1);
1653  GEOSGeom_destroy(g1);
1654 
1655  if (result == 2)
1656  {
1657  elog(ERROR,"GEOS isvalid() threw an error!");
1658  PG_RETURN_NULL(); /*never get here */
1659  }
1660 
1661  PG_FREE_IF_COPY(geom1, 0);
1662  PG_RETURN_BOOL(result);
1663 }
1664 
1665 PG_FUNCTION_INFO_V1(isvalidreason);
1666 Datum isvalidreason(PG_FUNCTION_ARGS)
1667 {
1668  GSERIALIZED *geom = NULL;
1669  char *reason_str = NULL;
1670  text *result = NULL;
1671  const GEOSGeometry *g1 = NULL;
1672 
1673  geom = PG_GETARG_GSERIALIZED_P(0);
1674 
1675  initGEOS(lwpgnotice, lwgeom_geos_error);
1676 
1677  g1 = POSTGIS2GEOS(geom);
1678  if ( g1 )
1679  {
1680  reason_str = GEOSisValidReason(g1);
1681  GEOSGeom_destroy((GEOSGeometry *)g1);
1682  if (!reason_str) HANDLE_GEOS_ERROR("GEOSisValidReason");
1683  result = cstring_to_text(reason_str);
1684  GEOSFree(reason_str);
1685  }
1686  else
1687  {
1688  result = cstring_to_text(lwgeom_geos_errmsg);
1689  }
1690 
1691  PG_FREE_IF_COPY(geom, 0);
1692  PG_RETURN_POINTER(result);
1693 }
1694 
1695 PG_FUNCTION_INFO_V1(isvaliddetail);
1696 Datum isvaliddetail(PG_FUNCTION_ARGS)
1697 {
1698  GSERIALIZED *geom = NULL;
1699  const GEOSGeometry *g1 = NULL;
1700  char *values[3]; /* valid bool, reason text, location geometry */
1701  char *geos_reason = NULL;
1702  char *reason = NULL;
1703  GEOSGeometry *geos_location = NULL;
1704  LWGEOM *location = NULL;
1705  char valid = 0;
1706  HeapTupleHeader result;
1707  TupleDesc tupdesc;
1708  HeapTuple tuple;
1709  AttInMetadata *attinmeta;
1710  int flags = 0;
1711 
1712  /*
1713  * Build a tuple description for a
1714  * valid_detail tuple
1715  */
1716  get_call_result_type(fcinfo, 0, &tupdesc);
1717  BlessTupleDesc(tupdesc);
1718 
1719  /*
1720  * generate attribute metadata needed later to produce
1721  * tuples from raw C strings
1722  */
1723  attinmeta = TupleDescGetAttInMetadata(tupdesc);
1724 
1725  geom = PG_GETARG_GSERIALIZED_P(0);
1726  flags = PG_GETARG_INT32(1);
1727 
1728  initGEOS(lwpgnotice, lwgeom_geos_error);
1729 
1730  g1 = POSTGIS2GEOS(geom);
1731 
1732  if ( g1 )
1733  {
1734  valid = GEOSisValidDetail(g1, flags, &geos_reason, &geos_location);
1735  GEOSGeom_destroy((GEOSGeometry *)g1);
1736  if ( geos_reason )
1737  {
1738  reason = pstrdup(geos_reason);
1739  GEOSFree(geos_reason);
1740  }
1741  if ( geos_location )
1742  {
1743  location = GEOS2LWGEOM(geos_location, GEOSHasZ(geos_location));
1744  GEOSGeom_destroy(geos_location);
1745  }
1746 
1747  if (valid == 2)
1748  {
1749  /* NOTE: should only happen on OOM or similar */
1750  lwpgerror("GEOS isvaliddetail() threw an exception!");
1751  PG_RETURN_NULL(); /* never gets here */
1752  }
1753  }
1754  else
1755  {
1756  /* TODO: check lwgeom_geos_errmsg for validity error */
1757  reason = pstrdup(lwgeom_geos_errmsg);
1758  }
1759 
1760  /* the boolean validity */
1761  values[0] = valid ? "t" : "f";
1762 
1763  /* the reason */
1764  values[1] = reason;
1765 
1766  /* the location */
1767  values[2] = location ? lwgeom_to_hexwkb_buffer(location, WKB_EXTENDED) : 0;
1768 
1769  tuple = BuildTupleFromCStrings(attinmeta, values);
1770  result = (HeapTupleHeader) palloc(tuple->t_len);
1771  memcpy(result, tuple->t_data, tuple->t_len);
1772  heap_freetuple(tuple);
1773 
1774  PG_RETURN_HEAPTUPLEHEADER(result);
1775 }
1776 
1784 PG_FUNCTION_INFO_V1(overlaps);
1785 Datum overlaps(PG_FUNCTION_ARGS)
1786 {
1787  GSERIALIZED *geom1;
1788  GSERIALIZED *geom2;
1789  GEOSGeometry *g1, *g2;
1790  char result;
1791  GBOX box1, box2;
1792 
1793  geom1 = PG_GETARG_GSERIALIZED_P(0);
1794  geom2 = PG_GETARG_GSERIALIZED_P(1);
1795  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1796 
1797  /* A.Overlaps(Empty) == FALSE */
1798  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1799  PG_RETURN_BOOL(false);
1800 
1801  /*
1802  * short-circuit 1: if geom2 bounding box does not overlap
1803  * geom1 bounding box we can return FALSE.
1804  */
1805  if ( gserialized_get_gbox_p(geom1, &box1) &&
1806  gserialized_get_gbox_p(geom2, &box2) )
1807  {
1808  if ( ! gbox_overlaps_2d(&box1, &box2) )
1809  {
1810  PG_RETURN_BOOL(false);
1811  }
1812  }
1813 
1814  initGEOS(lwpgnotice, lwgeom_geos_error);
1815 
1816  g1 = POSTGIS2GEOS(geom1);
1817  if (!g1)
1818  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1819 
1820  g2 = POSTGIS2GEOS(geom2);
1821 
1822  if (!g2)
1823  {
1824  GEOSGeom_destroy(g1);
1825  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1826  }
1827 
1828  result = GEOSOverlaps(g1,g2);
1829 
1830  GEOSGeom_destroy(g1);
1831  GEOSGeom_destroy(g2);
1832  if (result == 2) HANDLE_GEOS_ERROR("GEOSOverlaps");
1833 
1834  PG_FREE_IF_COPY(geom1, 0);
1835  PG_FREE_IF_COPY(geom2, 1);
1836 
1837  PG_RETURN_BOOL(result);
1838 }
1839 
1840 
1841 PG_FUNCTION_INFO_V1(contains);
1842 Datum contains(PG_FUNCTION_ARGS)
1843 {
1844  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
1845  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
1846  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
1847  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
1848  int result;
1849  GEOSGeometry *g1, *g2;
1850  GBOX box1, box2;
1851  PrepGeomCache *prep_cache;
1852  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1853 
1854  /* A.Contains(Empty) == FALSE */
1855  if (gserialized_is_empty(geom1) || gserialized_is_empty(geom2))
1856  PG_RETURN_BOOL(false);
1857 
1858  POSTGIS_DEBUG(3, "contains called.");
1859 
1860  /*
1861  ** short-circuit 1: if geom2 bounding box is not completely inside
1862  ** geom1 bounding box we can return FALSE.
1863  */
1864  if (gserialized_get_gbox_p(geom1, &box1) &&
1865  gserialized_get_gbox_p(geom2, &box2))
1866  {
1867  if (!gbox_contains_2d(&box1, &box2))
1868  PG_RETURN_BOOL(false);
1869  }
1870 
1871  /*
1872  ** short-circuit 2: if geom2 is a point and geom1 is a polygon
1873  ** call the point-in-polygon function.
1874  */
1875  if (is_poly(geom1) && is_point(geom2))
1876  {
1877  SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
1878  SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
1879  const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
1880  const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
1881  RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
1882  int retval;
1883 
1884  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1885  if (gserialized_get_type(gpoint) == POINTTYPE)
1886  {
1887  LWGEOM* point = lwgeom_from_gserialized(gpoint);
1888  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1889  lwgeom_free(point);
1890 
1891  retval = (pip_result == 1); /* completely inside */
1892  }
1893  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
1894  {
1895  LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
1896  int found_completely_inside = LW_FALSE;
1897 
1898  retval = LW_TRUE;
1899  for (uint32_t i = 0; i < mpoint->ngeoms; i++)
1900  {
1901  int pip_result;
1902  LWPOINT* pt = mpoint->geoms[i];
1903  /* We need to find at least one point that's completely inside the
1904  * polygons (pip_result == 1). As long as we have one point that's
1905  * completely inside, we can have as many as we want on the boundary
1906  * itself. (pip_result == 0)
1907  */
1908  if (lwpoint_is_empty(pt)) continue;
1909  pip_result = pip_short_circuit(cache, pt, gpoly);
1910  if (pip_result == 1)
1911  found_completely_inside = LW_TRUE;
1912 
1913  if (pip_result == -1) /* completely outside */
1914  {
1915  retval = LW_FALSE;
1916  break;
1917  }
1918  }
1919 
1920  retval = retval && found_completely_inside;
1921  lwmpoint_free(mpoint);
1922  }
1923  else
1924  {
1925  /* Never get here */
1926  elog(ERROR,"Type isn't point or multipoint!");
1927  PG_RETURN_BOOL(false);
1928  }
1929 
1930  return retval > 0;
1931  }
1932  else
1933  {
1934  POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
1935  }
1936 
1937  initGEOS(lwpgnotice, lwgeom_geos_error);
1938 
1939  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, NULL);
1940 
1941  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
1942  {
1943  g1 = POSTGIS2GEOS(geom2);
1944  if (!g1)
1945  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
1946 
1947  POSTGIS_DEBUG(4, "containsPrepared: cache is live, running preparedcontains");
1948  result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
1949  GEOSGeom_destroy(g1);
1950  }
1951  else
1952  {
1953  g1 = POSTGIS2GEOS(geom1);
1954  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1955  g2 = POSTGIS2GEOS(geom2);
1956  if (!g2)
1957  {
1958  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1959  GEOSGeom_destroy(g1);
1960  }
1961  POSTGIS_DEBUG(4, "containsPrepared: cache is not ready, running standard contains");
1962  result = GEOSContains( g1, g2);
1963  GEOSGeom_destroy(g1);
1964  GEOSGeom_destroy(g2);
1965  }
1966 
1967  if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1968 
1969  PG_RETURN_BOOL(result > 0);
1970 }
1971 
1972 
1973 PG_FUNCTION_INFO_V1(within);
1974 Datum within(PG_FUNCTION_ARGS)
1975 {
1976  PG_RETURN_DATUM(CallerFInfoFunctionCall2(contains, fcinfo->flinfo, InvalidOid,
1977  PG_GETARG_DATUM(1), PG_GETARG_DATUM(0)));
1978 }
1979 
1980 
1981 
1982 PG_FUNCTION_INFO_V1(containsproperly);
1983 Datum containsproperly(PG_FUNCTION_ARGS)
1984 {
1985  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
1986  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
1987  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
1988  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
1989  char result;
1990  GBOX box1, box2;
1991  PrepGeomCache * prep_cache;
1992 
1993  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
1994 
1995  /* A.ContainsProperly(Empty) == FALSE */
1996  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1997  PG_RETURN_BOOL(false);
1998 
1999  /*
2000  * short-circuit: if geom2 bounding box is not completely inside
2001  * geom1 bounding box we can return FALSE.
2002  */
2003  if ( gserialized_get_gbox_p(geom1, &box1) &&
2004  gserialized_get_gbox_p(geom2, &box2) )
2005  {
2006  if ( ! gbox_contains_2d(&box1, &box2) )
2007  PG_RETURN_BOOL(false);
2008  }
2009 
2010  initGEOS(lwpgnotice, lwgeom_geos_error);
2011 
2012  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
2013 
2014  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
2015  {
2016  GEOSGeometry *g = POSTGIS2GEOS(geom2);
2017  if (!g) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2018  result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
2019  GEOSGeom_destroy(g);
2020  }
2021  else
2022  {
2023  GEOSGeometry *g2;
2024  GEOSGeometry *g1;
2025 
2026  g1 = POSTGIS2GEOS(geom1);
2027  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2028  g2 = POSTGIS2GEOS(geom2);
2029  if (!g2)
2030  {
2031  GEOSGeom_destroy(g1);
2032  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2033  }
2034  result = GEOSRelatePattern( g1, g2, "T**FF*FF*" );
2035  GEOSGeom_destroy(g1);
2036  GEOSGeom_destroy(g2);
2037  }
2038 
2039  if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
2040 
2041  PG_RETURN_BOOL(result);
2042 }
2043 
2044 /*
2045  * Described at
2046  * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2047  */
2048 PG_FUNCTION_INFO_V1(covers);
2049 Datum covers(PG_FUNCTION_ARGS)
2050 {
2051  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2052  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2053  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2054  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2055  int result;
2056  GBOX box1, box2;
2057  PrepGeomCache *prep_cache;
2058 
2059 
2060  /* A.Covers(Empty) == FALSE */
2061  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2062  PG_RETURN_BOOL(false);
2063 
2064  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2065 
2066  /*
2067  * short-circuit 1: if geom2 bounding box is not completely inside
2068  * geom1 bounding box we can return FALSE.
2069  */
2070  if ( gserialized_get_gbox_p(geom1, &box1) &&
2071  gserialized_get_gbox_p(geom2, &box2) )
2072  {
2073  if ( ! gbox_contains_2d(&box1, &box2) )
2074  {
2075  PG_RETURN_BOOL(false);
2076  }
2077  }
2078  /*
2079  * short-circuit 2: if geom2 is a point and geom1 is a polygon
2080  * call the point-in-polygon function.
2081  */
2082  if (is_poly(geom1) && is_point(geom2))
2083  {
2084  SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2085  SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2086  const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2087  const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2088  RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2089  int retval;
2090 
2091  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2092  if (gserialized_get_type(gpoint) == POINTTYPE)
2093  {
2094  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2095  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2096  lwgeom_free(point);
2097 
2098  retval = (pip_result != -1); /* not outside */
2099  }
2100  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2101  {
2102  LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2103  uint32_t i;
2104 
2105  retval = LW_TRUE;
2106  for (i = 0; i < mpoint->ngeoms; i++)
2107  {
2108  LWPOINT *pt = mpoint->geoms[i];
2109  if (lwpoint_is_empty(pt)) continue;
2110  if (pip_short_circuit(cache, pt, gpoly) == -1)
2111  {
2112  retval = LW_FALSE;
2113  break;
2114  }
2115  }
2116 
2117  lwmpoint_free(mpoint);
2118  }
2119  else
2120  {
2121  /* Never get here */
2122  elog(ERROR,"Type isn't point or multipoint!");
2123  PG_RETURN_NULL();
2124  }
2125 
2126  PG_RETURN_BOOL(retval);
2127  }
2128  else
2129  {
2130  POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2131  }
2132 
2133  initGEOS(lwpgnotice, lwgeom_geos_error);
2134 
2135  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, 0);
2136 
2137  if ( prep_cache && prep_cache->prepared_geom && prep_cache->gcache.argnum == 1 )
2138  {
2139  GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
2140  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2141  result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
2142  GEOSGeom_destroy(g1);
2143  }
2144  else
2145  {
2146  GEOSGeometry *g1;
2147  GEOSGeometry *g2;
2148 
2149  g1 = POSTGIS2GEOS(geom1);
2150  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2151  g2 = POSTGIS2GEOS(geom2);
2152  if (!g2)
2153  {
2154  GEOSGeom_destroy(g1);
2155  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2156  }
2157  result = GEOSRelatePattern( g1, g2, "******FF*" );
2158  GEOSGeom_destroy(g1);
2159  GEOSGeom_destroy(g2);
2160  }
2161 
2162  if (result == 2) HANDLE_GEOS_ERROR("GEOSCovers");
2163 
2164  PG_RETURN_BOOL(result);
2165 }
2166 
2167 
2175 /*
2176  * Described at:
2177  * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2178  */
2179 PG_FUNCTION_INFO_V1(coveredby);
2180 Datum coveredby(PG_FUNCTION_ARGS)
2181 {
2182  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2183  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2184  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2185  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2186  GEOSGeometry *g1, *g2;
2187  int result;
2188  GBOX box1, box2;
2189  char *patt = "**F**F***";
2190 
2191  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2192 
2193  /* A.CoveredBy(Empty) == FALSE */
2194  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2195  PG_RETURN_BOOL(false);
2196 
2197  /*
2198  * short-circuit 1: if geom1 bounding box is not completely inside
2199  * geom2 bounding box we can return FALSE.
2200  */
2201  if ( gserialized_get_gbox_p(geom1, &box1) &&
2202  gserialized_get_gbox_p(geom2, &box2) )
2203  {
2204  if ( ! gbox_contains_2d(&box2, &box1) )
2205  {
2206  PG_RETURN_BOOL(false);
2207  }
2208 
2209  POSTGIS_DEBUG(3, "bounding box short-circuit missed.");
2210  }
2211  /*
2212  * short-circuit 2: if geom1 is a point and geom2 is a polygon
2213  * call the point-in-polygon function.
2214  */
2215  if (is_point(geom1) && is_poly(geom2))
2216  {
2217  SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2218  SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2219  const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2220  const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2221  RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2222  int retval;
2223 
2224  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2225  if (gserialized_get_type(gpoint) == POINTTYPE)
2226  {
2227  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2228  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2229  lwgeom_free(point);
2230 
2231  retval = (pip_result != -1); /* not outside */
2232  }
2233  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2234  {
2235  LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2236  uint32_t i;
2237 
2238  retval = LW_TRUE;
2239  for (i = 0; i < mpoint->ngeoms; i++)
2240  {
2241  LWPOINT *pt = mpoint->geoms[i];
2242  if (lwpoint_is_empty(pt)) continue;
2243  if (pip_short_circuit(cache, pt, gpoly) == -1)
2244  {
2245  retval = LW_FALSE;
2246  break;
2247  }
2248  }
2249 
2250  lwmpoint_free(mpoint);
2251  }
2252  else
2253  {
2254  /* Never get here */
2255  elog(ERROR,"Type isn't point or multipoint!");
2256  PG_RETURN_NULL();
2257  }
2258 
2259  PG_RETURN_BOOL(retval);
2260  }
2261  else
2262  {
2263  POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2264  }
2265 
2266  initGEOS(lwpgnotice, lwgeom_geos_error);
2267 
2268  g1 = POSTGIS2GEOS(geom1);
2269 
2270  if (!g1)
2271  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2272 
2273  g2 = POSTGIS2GEOS(geom2);
2274 
2275  if (!g2)
2276  {
2277  GEOSGeom_destroy(g1);
2278  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2279  }
2280 
2281  result = GEOSRelatePattern(g1,g2,patt);
2282 
2283  GEOSGeom_destroy(g1);
2284  GEOSGeom_destroy(g2);
2285 
2286  if (result == 2) HANDLE_GEOS_ERROR("GEOSCoveredBy");
2287 
2288  PG_RETURN_BOOL(result);
2289 }
2290 
2291 PG_FUNCTION_INFO_V1(crosses);
2292 Datum crosses(PG_FUNCTION_ARGS)
2293 {
2294  GSERIALIZED *geom1;
2295  GSERIALIZED *geom2;
2296  GEOSGeometry *g1, *g2;
2297  int result;
2298  GBOX box1, box2;
2299 
2300  geom1 = PG_GETARG_GSERIALIZED_P(0);
2301  geom2 = PG_GETARG_GSERIALIZED_P(1);
2302  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2303 
2304  /* A.Crosses(Empty) == FALSE */
2305  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2306  PG_RETURN_BOOL(false);
2307 
2308  /*
2309  * short-circuit 1: if geom2 bounding box does not overlap
2310  * geom1 bounding box we can return FALSE.
2311  */
2312  if ( gserialized_get_gbox_p(geom1, &box1) &&
2313  gserialized_get_gbox_p(geom2, &box2) )
2314  {
2315  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2316  {
2317  PG_RETURN_BOOL(false);
2318  }
2319  }
2320 
2321  initGEOS(lwpgnotice, lwgeom_geos_error);
2322 
2323  g1 = POSTGIS2GEOS(geom1);
2324  if (!g1)
2325  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2326 
2327  g2 = POSTGIS2GEOS(geom2);
2328  if (!g2)
2329  {
2330  GEOSGeom_destroy(g1);
2331  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2332  }
2333 
2334  result = GEOSCrosses(g1,g2);
2335 
2336  GEOSGeom_destroy(g1);
2337  GEOSGeom_destroy(g2);
2338 
2339  if (result == 2) HANDLE_GEOS_ERROR("GEOSCrosses");
2340 
2341  PG_FREE_IF_COPY(geom1, 0);
2342  PG_FREE_IF_COPY(geom2, 1);
2343 
2344  PG_RETURN_BOOL(result);
2345 }
2346 
2347 PG_FUNCTION_INFO_V1(ST_Intersects);
2348 Datum ST_Intersects(PG_FUNCTION_ARGS)
2349 {
2350  SHARED_GSERIALIZED *shared_geom1 = ToastCacheGetGeometry(fcinfo, 0);
2351  SHARED_GSERIALIZED *shared_geom2 = ToastCacheGetGeometry(fcinfo, 1);
2352  const GSERIALIZED *geom1 = shared_gserialized_get(shared_geom1);
2353  const GSERIALIZED *geom2 = shared_gserialized_get(shared_geom2);
2354  int result;
2355  GBOX box1, box2;
2356  PrepGeomCache *prep_cache;
2357 
2358  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2359 
2360  /* A.Intersects(Empty) == FALSE */
2361  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2362  PG_RETURN_BOOL(false);
2363 
2364  /*
2365  * short-circuit 1: if geom2 bounding box does not overlap
2366  * geom1 bounding box we can return FALSE.
2367  */
2368  if ( gserialized_get_gbox_p(geom1, &box1) &&
2369  gserialized_get_gbox_p(geom2, &box2) )
2370  {
2371  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2372  PG_RETURN_BOOL(false);
2373  }
2374 
2375  /*
2376  * short-circuit 2: if the geoms are a point and a polygon,
2377  * call the point_outside_polygon function.
2378  */
2379  if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
2380  {
2381  SHARED_GSERIALIZED *shared_gpoly = is_poly(geom1) ? shared_geom1 : shared_geom2;
2382  SHARED_GSERIALIZED *shared_gpoint = is_point(geom1) ? shared_geom1 : shared_geom2;
2383  const GSERIALIZED *gpoly = shared_gserialized_get(shared_gpoly);
2384  const GSERIALIZED *gpoint = shared_gserialized_get(shared_gpoint);
2385  RTREE_POLY_CACHE *cache = GetRtreeCache(fcinfo, shared_gpoly);
2386  int retval;
2387 
2388  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2389  if (gserialized_get_type(gpoint) == POINTTYPE)
2390  {
2391  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2392  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2393  lwgeom_free(point);
2394 
2395  retval = (pip_result != -1); /* not outside */
2396  }
2397  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2398  {
2399  LWMPOINT* mpoint = lwgeom_as_lwmpoint(lwgeom_from_gserialized(gpoint));
2400  retval = LW_FALSE;
2401  for (uint32_t i = 0; i < mpoint->ngeoms; i++)
2402  {
2403  LWPOINT *pt = mpoint->geoms[i];
2404  if (lwpoint_is_empty(pt)) continue;
2405  if (pip_short_circuit(cache, pt, gpoly) != -1) /* not outside */
2406  {
2407  retval = LW_TRUE;
2408  break;
2409  }
2410  }
2411 
2412  lwmpoint_free(mpoint);
2413  }
2414  else
2415  {
2416  /* Never get here */
2417  elog(ERROR,"Type isn't point or multipoint!");
2418  PG_RETURN_NULL();
2419  }
2420 
2421  PG_RETURN_BOOL(retval);
2422  }
2423 
2424  initGEOS(lwpgnotice, lwgeom_geos_error);
2425  prep_cache = GetPrepGeomCache(fcinfo, shared_geom1, shared_geom2);
2426 
2427  if ( prep_cache && prep_cache->prepared_geom )
2428  {
2429  if ( prep_cache->gcache.argnum == 1 )
2430  {
2431  GEOSGeometry *g = POSTGIS2GEOS(geom2);
2432  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2433  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2434  GEOSGeom_destroy(g);
2435  }
2436  else
2437  {
2438  GEOSGeometry *g = POSTGIS2GEOS(geom1);
2439  if (!g)
2440  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2441  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2442  GEOSGeom_destroy(g);
2443  }
2444  }
2445  else
2446  {
2447  GEOSGeometry *g1;
2448  GEOSGeometry *g2;
2449  g1 = POSTGIS2GEOS(geom1);
2450  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2451  g2 = POSTGIS2GEOS(geom2);
2452  if (!g2)
2453  {
2454  GEOSGeom_destroy(g1);
2455  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2456  }
2457  result = GEOSIntersects( g1, g2);
2458  GEOSGeom_destroy(g1);
2459  GEOSGeom_destroy(g2);
2460  }
2461 
2462  if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
2463 
2464  PG_RETURN_BOOL(result);
2465 }
2466 
2467 
2468 PG_FUNCTION_INFO_V1(touches);
2469 Datum touches(PG_FUNCTION_ARGS)
2470 {
2471  GSERIALIZED *geom1;
2472  GSERIALIZED *geom2;
2473  GEOSGeometry *g1, *g2;
2474  char result;
2475  GBOX box1, box2;
2476 
2477  geom1 = PG_GETARG_GSERIALIZED_P(0);
2478  geom2 = PG_GETARG_GSERIALIZED_P(1);
2479  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2480 
2481  /* A.Touches(Empty) == FALSE */
2482  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2483  PG_RETURN_BOOL(false);
2484 
2485  /*
2486  * short-circuit 1: if geom2 bounding box does not overlap
2487  * geom1 bounding box we can return FALSE.
2488  */
2489  if ( gserialized_get_gbox_p(geom1, &box1) &&
2490  gserialized_get_gbox_p(geom2, &box2) )
2491  {
2492  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2493  {
2494  PG_RETURN_BOOL(false);
2495  }
2496  }
2497 
2498  initGEOS(lwpgnotice, lwgeom_geos_error);
2499 
2500  g1 = POSTGIS2GEOS(geom1 );
2501  if (!g1)
2502  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2503 
2504  g2 = POSTGIS2GEOS(geom2 );
2505  if (!g2)
2506  {
2507  GEOSGeom_destroy(g1);
2508  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2509  }
2510 
2511  result = GEOSTouches(g1,g2);
2512 
2513  GEOSGeom_destroy(g1);
2514  GEOSGeom_destroy(g2);
2515 
2516  if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
2517 
2518  PG_FREE_IF_COPY(geom1, 0);
2519  PG_FREE_IF_COPY(geom2, 1);
2520 
2521  PG_RETURN_BOOL(result);
2522 }
2523 
2524 
2525 PG_FUNCTION_INFO_V1(disjoint);
2526 Datum disjoint(PG_FUNCTION_ARGS)
2527 {
2528  GSERIALIZED *geom1;
2529  GSERIALIZED *geom2;
2530  GEOSGeometry *g1, *g2;
2531  char result;
2532  GBOX box1, box2;
2533 
2534  geom1 = PG_GETARG_GSERIALIZED_P(0);
2535  geom2 = PG_GETARG_GSERIALIZED_P(1);
2536  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2537 
2538  /* A.Disjoint(Empty) == TRUE */
2539  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2540  PG_RETURN_BOOL(true);
2541 
2542  /*
2543  * short-circuit 1: if geom2 bounding box does not overlap
2544  * geom1 bounding box we can return TRUE.
2545  */
2546  if ( gserialized_get_gbox_p(geom1, &box1) &&
2547  gserialized_get_gbox_p(geom2, &box2) )
2548  {
2549  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2550  {
2551  PG_RETURN_BOOL(true);
2552  }
2553  }
2554 
2555  initGEOS(lwpgnotice, lwgeom_geos_error);
2556 
2557  g1 = POSTGIS2GEOS(geom1);
2558  if (!g1)
2559  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2560 
2561  g2 = POSTGIS2GEOS(geom2);
2562  if (!g2)
2563  {
2564  GEOSGeom_destroy(g1);
2565  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2566  }
2567 
2568  result = GEOSDisjoint(g1,g2);
2569 
2570  GEOSGeom_destroy(g1);
2571  GEOSGeom_destroy(g2);
2572 
2573  if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
2574 
2575  PG_FREE_IF_COPY(geom1, 0);
2576  PG_FREE_IF_COPY(geom2, 1);
2577 
2578  PG_RETURN_BOOL(result);
2579 }
2580 
2581 
2582 PG_FUNCTION_INFO_V1(relate_pattern);
2583 Datum relate_pattern(PG_FUNCTION_ARGS)
2584 {
2585  GSERIALIZED *geom1;
2586  GSERIALIZED *geom2;
2587  char *patt;
2588  text *ptxt;
2589  char result;
2590  GEOSGeometry *g1, *g2;
2591  size_t i;
2592 
2593  geom1 = PG_GETARG_GSERIALIZED_P(0);
2594  geom2 = PG_GETARG_GSERIALIZED_P(1);
2595  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2596 
2597  /* TODO handle empty */
2598 
2599  initGEOS(lwpgnotice, lwgeom_geos_error);
2600 
2601  g1 = POSTGIS2GEOS(geom1);
2602  if (!g1)
2603  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2604  g2 = POSTGIS2GEOS(geom2);
2605  if (!g2)
2606  {
2607  GEOSGeom_destroy(g1);
2608  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2609  }
2610 
2611  /* Ensure DE9IM pattern is no more than 9 chars */
2612  ptxt = DatumGetTextP(DirectFunctionCall2(text_left,
2613  PG_GETARG_DATUM(2), Int32GetDatum(9)));
2614  patt = text_to_cstring(ptxt);
2615 
2616  /*
2617  ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2618  */
2619  for ( i = 0; i < strlen(patt); i++ )
2620  {
2621  if ( patt[i] == 't' ) patt[i] = 'T';
2622  if ( patt[i] == 'f' ) patt[i] = 'F';
2623  }
2624 
2625  result = GEOSRelatePattern(g1,g2,patt);
2626  GEOSGeom_destroy(g1);
2627  GEOSGeom_destroy(g2);
2628  pfree(patt);
2629 
2630  if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
2631 
2632  PG_FREE_IF_COPY(geom1, 0);
2633  PG_FREE_IF_COPY(geom2, 1);
2634 
2635  PG_RETURN_BOOL(result);
2636 }
2637 
2638 
2639 
2640 PG_FUNCTION_INFO_V1(relate_full);
2641 Datum relate_full(PG_FUNCTION_ARGS)
2642 {
2643  GSERIALIZED *geom1;
2644  GSERIALIZED *geom2;
2645  GEOSGeometry *g1, *g2;
2646  char *relate_str;
2647  text *result;
2648  int bnr = GEOSRELATE_BNR_OGC;
2649 
2650  /* TODO handle empty */
2651  geom1 = PG_GETARG_GSERIALIZED_P(0);
2652  geom2 = PG_GETARG_GSERIALIZED_P(1);
2653  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2654 
2655  if ( PG_NARGS() > 2 )
2656  bnr = PG_GETARG_INT32(2);
2657 
2658  initGEOS(lwpgnotice, lwgeom_geos_error);
2659 
2660  g1 = POSTGIS2GEOS(geom1 );
2661  if (!g1)
2662  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2663  g2 = POSTGIS2GEOS(geom2 );
2664  if (!g2)
2665  {
2666  GEOSGeom_destroy(g1);
2667  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2668  }
2669 
2670  POSTGIS_DEBUG(3, "constructed geometries ");
2671 
2672  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2673  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2674 
2675  relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2676 
2677  GEOSGeom_destroy(g1);
2678  GEOSGeom_destroy(g2);
2679 
2680  if (!relate_str) HANDLE_GEOS_ERROR("GEOSRelate");
2681 
2682  result = cstring_to_text(relate_str);
2683  GEOSFree(relate_str);
2684 
2685  PG_FREE_IF_COPY(geom1, 0);
2686  PG_FREE_IF_COPY(geom2, 1);
2687 
2688  PG_RETURN_TEXT_P(result);
2689 }
2690 
2691 
2692 PG_FUNCTION_INFO_V1(ST_Equals);
2693 Datum ST_Equals(PG_FUNCTION_ARGS)
2694 {
2695  GSERIALIZED *geom1;
2696  GSERIALIZED *geom2;
2697  GEOSGeometry *g1, *g2;
2698  char result;
2699  GBOX box1, box2;
2700 
2701  geom1 = PG_GETARG_GSERIALIZED_P(0);
2702  geom2 = PG_GETARG_GSERIALIZED_P(1);
2703  gserialized_error_if_srid_mismatch(geom1, geom2, __func__);
2704 
2705  /* Empty == Empty */
2706  if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
2707  PG_RETURN_BOOL(true);
2708 
2709  /*
2710  * short-circuit: If geom1 and geom2 do not have the same bounding box
2711  * we can return FALSE.
2712  */
2713  if ( gserialized_get_gbox_p(geom1, &box1) &&
2714  gserialized_get_gbox_p(geom2, &box2) )
2715  {
2716  if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
2717  {
2718  PG_RETURN_BOOL(false);
2719  }
2720  }
2721 
2722  /*
2723  * short-circuit: if geom1 and geom2 are binary-equivalent, we can return
2724  * TRUE. This is much faster than doing the comparison using GEOS.
2725  */
2726  if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
2727  PG_RETURN_BOOL(true);
2728  }
2729 
2730  initGEOS(lwpgnotice, lwgeom_geos_error);
2731 
2732  g1 = POSTGIS2GEOS(geom1);
2733 
2734  if (!g1)
2735  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2736 
2737  g2 = POSTGIS2GEOS(geom2);
2738 
2739  if (!g2)
2740  {
2741  GEOSGeom_destroy(g1);
2742  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2743  }
2744 
2745  result = GEOSEquals(g1,g2);
2746 
2747  GEOSGeom_destroy(g1);
2748  GEOSGeom_destroy(g2);
2749 
2750  if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
2751 
2752  PG_FREE_IF_COPY(geom1, 0);
2753  PG_FREE_IF_COPY(geom2, 1);
2754 
2755  PG_RETURN_BOOL(result);
2756 }
2757 
2758 PG_FUNCTION_INFO_V1(issimple);
2759 Datum issimple(PG_FUNCTION_ARGS)
2760 {
2761  GSERIALIZED *geom;
2762  LWGEOM *lwgeom_in;
2763  int result;
2764 
2765  POSTGIS_DEBUG(2, "issimple called");
2766 
2767  geom = PG_GETARG_GSERIALIZED_P(0);
2768 
2769  if ( gserialized_is_empty(geom) )
2770  PG_RETURN_BOOL(true);
2771 
2772  lwgeom_in = lwgeom_from_gserialized(geom);
2773  result = lwgeom_is_simple(lwgeom_in);
2774  lwgeom_free(lwgeom_in) ;
2775  PG_FREE_IF_COPY(geom, 0);
2776 
2777  if (result == -1) {
2778  PG_RETURN_NULL(); /*never get here */
2779  }
2780 
2781  PG_RETURN_BOOL(result);
2782 }
2783 
2784 PG_FUNCTION_INFO_V1(isring);
2785 Datum isring(PG_FUNCTION_ARGS)
2786 {
2787  GSERIALIZED *geom;
2788  GEOSGeometry *g1;
2789  int result;
2790 
2791  geom = PG_GETARG_GSERIALIZED_P(0);
2792 
2793  /* Empty things can't close */
2794  if ( gserialized_is_empty(geom) )
2795  PG_RETURN_BOOL(false);
2796 
2797  initGEOS(lwpgnotice, lwgeom_geos_error);
2798 
2799  g1 = POSTGIS2GEOS(geom);
2800  if (!g1)
2801  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2802 
2803  if ( GEOSGeomTypeId(g1) != GEOS_LINESTRING )
2804  {
2805  GEOSGeom_destroy(g1);
2806  elog(ERROR, "ST_IsRing() should only be called on a linear feature");
2807  }
2808 
2809  result = GEOSisRing(g1);
2810  GEOSGeom_destroy(g1);
2811 
2812  if (result == 2) HANDLE_GEOS_ERROR("GEOSisRing");
2813 
2814  PG_FREE_IF_COPY(geom, 0);
2815  PG_RETURN_BOOL(result);
2816 }
2817 
2818 GSERIALIZED *
2819 GEOS2POSTGIS(GEOSGeom geom, char want3d)
2820 {
2821  LWGEOM *lwgeom;
2822  GSERIALIZED *result;
2823 
2824  lwgeom = GEOS2LWGEOM(geom, want3d);
2825  if ( ! lwgeom )
2826  {
2827  lwpgerror("%s: GEOS2LWGEOM returned NULL", __func__);
2828  return NULL;
2829  }
2830 
2831  POSTGIS_DEBUGF(4, "%s: GEOS2LWGEOM returned a %s", __func__, lwgeom_summary(lwgeom, 0));
2832 
2833  if (lwgeom_needs_bbox(lwgeom)) lwgeom_add_bbox(lwgeom);
2834 
2835  result = geometry_serialize(lwgeom);
2836  lwgeom_free(lwgeom);
2837 
2838  return result;
2839 }
2840 
2841 /*-----=POSTGIS2GEOS= */
2842 
2843 GEOSGeometry *
2844 POSTGIS2GEOS(const GSERIALIZED *pglwgeom)
2845 {
2846  GEOSGeometry *ret;
2847  LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
2848  if ( ! lwgeom )
2849  {
2850  lwpgerror("POSTGIS2GEOS: unable to deserialize input");
2851  return NULL;
2852  }
2853  ret = LWGEOM2GEOS(lwgeom, 0);
2854  lwgeom_free(lwgeom);
2855 
2856  return ret;
2857 }
2858 
2859 uint32_t array_nelems_not_null(ArrayType* array) {
2860  ArrayIterator iterator;
2861  Datum value;
2862  bool isnull;
2863  uint32_t nelems_not_null = 0;
2864  iterator = array_create_iterator(array, 0, NULL);
2865  while(array_iterate(iterator, &value, &isnull) )
2866  if (!isnull)
2867  nelems_not_null++;
2868 
2869  array_free_iterator(iterator);
2870 
2871  return nelems_not_null;
2872 }
2873 
2874 /* ARRAY2LWGEOM: Converts the non-null elements of a Postgres array into a LWGEOM* array */
2875 LWGEOM** ARRAY2LWGEOM(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2876 {
2877  ArrayIterator iterator;
2878  Datum value;
2879  bool isnull;
2880  bool gotsrid = false;
2881  uint32_t i = 0;
2882 
2883  LWGEOM** lw_geoms = palloc(nelems * sizeof(LWGEOM*));
2884 
2885  iterator = array_create_iterator(array, 0, NULL);
2886 
2887  while (array_iterate(iterator, &value, &isnull))
2888  {
2889  GSERIALIZED *geom = (GSERIALIZED *)DatumGetPointer(value);
2890 
2891  if (isnull)
2892  continue;
2893 
2894  *is3d = *is3d || gserialized_has_z(geom);
2895 
2896  lw_geoms[i] = lwgeom_from_gserialized(geom);
2897  if (!lw_geoms[i]) /* error in creation */
2898  {
2899  lwpgerror("Geometry deserializing geometry");
2900  return NULL;
2901  }
2902  if (!gotsrid)
2903  {
2904  gotsrid = true;
2905  *srid = gserialized_get_srid(geom);
2906  }
2907  else
2908  gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2909 
2910  i++;
2911  }
2912 
2913  return lw_geoms;
2914 }
2915 
2916 /* ARRAY2GEOS: Converts the non-null elements of a Postgres array into a GEOSGeometry* array */
2917 GEOSGeometry** ARRAY2GEOS(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2918 {
2919  ArrayIterator iterator;
2920  Datum value;
2921  bool isnull;
2922  bool gotsrid = false;
2923  uint32_t i = 0;
2924 
2925  GEOSGeometry** geos_geoms = palloc(nelems * sizeof(GEOSGeometry*));
2926 
2927  iterator = array_create_iterator(array, 0, NULL);
2928 
2929  while(array_iterate(iterator, &value, &isnull))
2930  {
2931  GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2932 
2933  if (isnull)
2934  continue;
2935 
2936  *is3d = *is3d || gserialized_has_z(geom);
2937 
2938  geos_geoms[i] = POSTGIS2GEOS(geom);
2939  if (!geos_geoms[i])
2940  {
2941  uint32_t j;
2942  lwpgerror("Geometry could not be converted to GEOS");
2943 
2944  for (j = 0; j < i; j++) {
2945  GEOSGeom_destroy(geos_geoms[j]);
2946  }
2947  return NULL;
2948  }
2949 
2950  if (!gotsrid)
2951  {
2952  *srid = gserialized_get_srid(geom);
2953  gotsrid = true;
2954  }
2955  else if (*srid != gserialized_get_srid(geom))
2956  {
2957  uint32_t j;
2958  for (j = 0; j <= i; j++) {
2959  GEOSGeom_destroy(geos_geoms[j]);
2960  }
2961  gserialized_error_if_srid_mismatch_reference(geom, *srid, __func__);
2962  return NULL;
2963  }
2964 
2965  i++;
2966  }
2967 
2968  array_free_iterator(iterator);
2969  return geos_geoms;
2970 }
2971 
2972 PG_FUNCTION_INFO_V1(GEOSnoop);
2973 Datum GEOSnoop(PG_FUNCTION_ARGS)
2974 {
2975  GSERIALIZED *geom;
2976  GEOSGeometry *geosgeom;
2977  GSERIALIZED *lwgeom_result;
2978 
2979  initGEOS(lwpgnotice, lwgeom_geos_error);
2980 
2981  geom = PG_GETARG_GSERIALIZED_P(0);
2982  geosgeom = POSTGIS2GEOS(geom);
2983  if ( ! geosgeom ) PG_RETURN_NULL();
2984 
2985  lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
2986  GEOSGeom_destroy(geosgeom);
2987 
2988  PG_FREE_IF_COPY(geom, 0);
2989 
2990  PG_RETURN_POINTER(lwgeom_result);
2991 }
2992 
2993 PG_FUNCTION_INFO_V1(polygonize_garray);
2994 Datum polygonize_garray(PG_FUNCTION_ARGS)
2995 {
2996  ArrayType *array;
2997  int is3d = 0;
2998  uint32 nelems, i;
2999  GSERIALIZED *result;
3000  GEOSGeometry *geos_result;
3001  const GEOSGeometry **vgeoms;
3002  int32_t srid = SRID_UNKNOWN;
3003 #if POSTGIS_DEBUG_LEVEL >= 3
3004  static int call=1;
3005 #endif
3006 
3007 #if POSTGIS_DEBUG_LEVEL >= 3
3008  call++;
3009 #endif
3010 
3011  if (PG_ARGISNULL(0))
3012  PG_RETURN_NULL();
3013 
3014  array = PG_GETARG_ARRAYTYPE_P(0);
3015  nelems = array_nelems_not_null(array);
3016 
3017  if (nelems == 0)
3018  PG_RETURN_NULL();
3019 
3020  POSTGIS_DEBUGF(3, "polygonize_garray: number of non-null elements: %d", nelems);
3021 
3022  /* Ok, we really need geos now ;) */
3023  initGEOS(lwpgnotice, lwgeom_geos_error);
3024 
3025  vgeoms = (const GEOSGeometry**) ARRAY2GEOS(array, nelems, &is3d, &srid);
3026 
3027  POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
3028 
3029  geos_result = GEOSPolygonize(vgeoms, nelems);
3030 
3031  POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
3032 
3033  for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
3034  pfree(vgeoms);
3035 
3036  if ( ! geos_result ) PG_RETURN_NULL();
3037 
3038  GEOSSetSRID(geos_result, srid);
3039  result = GEOS2POSTGIS(geos_result, is3d);
3040  GEOSGeom_destroy(geos_result);
3041  if (!result)
3042  {
3043  elog(ERROR, "%s returned an error", __func__);
3044  PG_RETURN_NULL(); /*never get here */
3045  }
3046 
3047  PG_RETURN_POINTER(result);
3048 }
3049 
3050 
3051 PG_FUNCTION_INFO_V1(clusterintersecting_garray);
3052 Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
3053 {
3054  Datum* result_array_data;
3055  ArrayType *array, *result;
3056  int is3d = 0;
3057  uint32 nelems, nclusters, i;
3058  GEOSGeometry **geos_inputs, **geos_results;
3059  int32_t srid = SRID_UNKNOWN;
3060 
3061  /* Parameters used to construct a result array */
3062  int16 elmlen;
3063  bool elmbyval;
3064  char elmalign;
3065 
3066  /* Null array, null geometry (should be empty?) */
3067  if (PG_ARGISNULL(0))
3068  PG_RETURN_NULL();
3069 
3070  array = PG_GETARG_ARRAYTYPE_P(0);
3071  nelems = array_nelems_not_null(array);
3072 
3073  POSTGIS_DEBUGF(3, "clusterintersecting_garray: number of non-null elements: %d", nelems);
3074 
3075  if ( nelems == 0 ) PG_RETURN_NULL();
3076 
3077  /* TODO short-circuit for one element? */
3078 
3079  /* Ok, we really need geos now ;) */
3080  initGEOS(lwpgnotice, lwgeom_geos_error);
3081 
3082  geos_inputs = ARRAY2GEOS(array, nelems, &is3d, &srid);
3083  if(!geos_inputs)
3084  {
3085  PG_RETURN_NULL();
3086  }
3087 
3088  if (cluster_intersecting(geos_inputs, nelems, &geos_results, &nclusters) != LW_SUCCESS)
3089  {
3090  elog(ERROR, "clusterintersecting: Error performing clustering");
3091  PG_RETURN_NULL();
3092  }
3093  pfree(geos_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3094 
3095  if (!geos_results) PG_RETURN_NULL();
3096 
3097  result_array_data = palloc(nclusters * sizeof(Datum));
3098  for (i=0; i<nclusters; ++i)
3099  {
3100  result_array_data[i] = PointerGetDatum(GEOS2POSTGIS(geos_results[i], is3d));
3101  GEOSGeom_destroy(geos_results[i]);
3102  }
3103  lwfree(geos_results);
3104 
3105  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3106  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3107 
3108  if (!result)
3109  {
3110  elog(ERROR, "clusterintersecting: Error constructing return-array");
3111  PG_RETURN_NULL();
3112  }
3113 
3114  PG_RETURN_POINTER(result);
3115 }
3116 
3117 PG_FUNCTION_INFO_V1(cluster_within_distance_garray);
3118 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS)
3119 {
3120  Datum* result_array_data;
3121  ArrayType *array, *result;
3122  int is3d = 0;
3123  uint32 nelems, nclusters, i;
3124  LWGEOM** lw_inputs;
3125  LWGEOM** lw_results;
3126  double tolerance;
3127  int32_t srid = SRID_UNKNOWN;
3128 
3129  /* Parameters used to construct a result array */
3130  int16 elmlen;
3131  bool elmbyval;
3132  char elmalign;
3133 
3134  /* Null array, null geometry (should be empty?) */
3135  if (PG_ARGISNULL(0))
3136  PG_RETURN_NULL();
3137 
3138  array = PG_GETARG_ARRAYTYPE_P(0);
3139 
3140  tolerance = PG_GETARG_FLOAT8(1);
3141  if (tolerance < 0)
3142  {
3143  lwpgerror("Tolerance must be a positive number.");
3144  PG_RETURN_NULL();
3145  }
3146 
3147  nelems = array_nelems_not_null(array);
3148 
3149  POSTGIS_DEBUGF(3, "cluster_within_distance_garray: number of non-null elements: %d", nelems);
3150 
3151  if ( nelems == 0 ) PG_RETURN_NULL();
3152 
3153  /* TODO short-circuit for one element? */
3154 
3155  /* Ok, we really need geos now ;) */
3156  initGEOS(lwpgnotice, lwgeom_geos_error);
3157 
3158  lw_inputs = ARRAY2LWGEOM(array, nelems, &is3d, &srid);
3159  if (!lw_inputs)
3160  {
3161  PG_RETURN_NULL();
3162  }
3163 
3164  if (cluster_within_distance(lw_inputs, nelems, tolerance, &lw_results, &nclusters) != LW_SUCCESS)
3165  {
3166  elog(ERROR, "cluster_within: Error performing clustering");
3167  PG_RETURN_NULL();
3168  }
3169  pfree(lw_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3170 
3171  if (!lw_results) PG_RETURN_NULL();
3172 
3173  result_array_data = palloc(nclusters * sizeof(Datum));
3174  for (i=0; i<nclusters; ++i)
3175  {
3176  result_array_data[i] = PointerGetDatum(geometry_serialize(lw_results[i]));
3177  lwgeom_free(lw_results[i]);
3178  }
3179  lwfree(lw_results);
3180 
3181  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3182  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3183 
3184  if (!result)
3185  {
3186  elog(ERROR, "clusterwithin: Error constructing return-array");
3187  PG_RETURN_NULL();
3188  }
3189 
3190  PG_RETURN_POINTER(result);
3191 }
3192 
3193 PG_FUNCTION_INFO_V1(linemerge);
3194 Datum linemerge(PG_FUNCTION_ARGS)
3195 {
3196  GSERIALIZED *geom1;
3197  GSERIALIZED *result;
3198  LWGEOM *lwgeom1, *lwresult ;
3199 
3200  geom1 = PG_GETARG_GSERIALIZED_P(0);
3201 
3202 
3203  lwgeom1 = lwgeom_from_gserialized(geom1) ;
3204 
3205  lwresult = lwgeom_linemerge(lwgeom1);
3206  result = geometry_serialize(lwresult) ;
3207 
3208  lwgeom_free(lwgeom1) ;
3209  lwgeom_free(lwresult) ;
3210 
3211  PG_FREE_IF_COPY(geom1, 0);
3212 
3213  PG_RETURN_POINTER(result);
3214 }
3215 
3216 /*
3217  * Take a geometry and return an areal geometry
3218  * (Polygon or MultiPolygon).
3219  * Actually a wrapper around GEOSpolygonize,
3220  * transforming the resulting collection into
3221  * a valid polygon Geometry.
3222  */
3223 PG_FUNCTION_INFO_V1(ST_BuildArea);
3224 Datum ST_BuildArea(PG_FUNCTION_ARGS)
3225 {
3226  GSERIALIZED *result;
3227  GSERIALIZED *geom;
3228  LWGEOM *lwgeom_in, *lwgeom_out;
3229 
3230  geom = PG_GETARG_GSERIALIZED_P(0);
3231  lwgeom_in = lwgeom_from_gserialized(geom);
3232 
3233  lwgeom_out = lwgeom_buildarea(lwgeom_in);
3234  lwgeom_free(lwgeom_in) ;
3235 
3236  if ( ! lwgeom_out )
3237  {
3238  PG_FREE_IF_COPY(geom, 0);
3239  PG_RETURN_NULL();
3240  }
3241 
3242  result = geometry_serialize(lwgeom_out) ;
3243  lwgeom_free(lwgeom_out) ;
3244 
3245  PG_FREE_IF_COPY(geom, 0);
3246  PG_RETURN_POINTER(result);
3247 }
3248 
3249 /*
3250  * Take the vertices of a geometry and builds
3251  * Delaunay triangles around them.
3252  */
3253 PG_FUNCTION_INFO_V1(ST_DelaunayTriangles);
3254 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
3255 {
3256  GSERIALIZED *result;
3257  GSERIALIZED *geom;
3258  LWGEOM *lwgeom_in, *lwgeom_out;
3259  double tolerance = 0.0;
3260  int flags = 0;
3261 
3262  geom = PG_GETARG_GSERIALIZED_P(0);
3263  tolerance = PG_GETARG_FLOAT8(1);
3264  flags = PG_GETARG_INT32(2);
3265 
3266  lwgeom_in = lwgeom_from_gserialized(geom);
3267  lwgeom_out = lwgeom_delaunay_triangulation(lwgeom_in, tolerance, flags);
3268  lwgeom_free(lwgeom_in) ;
3269 
3270  if ( ! lwgeom_out )
3271  {
3272  PG_FREE_IF_COPY(geom, 0);
3273  PG_RETURN_NULL();
3274  }
3275 
3276  result = geometry_serialize(lwgeom_out) ;
3277  lwgeom_free(lwgeom_out) ;
3278 
3279  PG_FREE_IF_COPY(geom, 0);
3280  PG_RETURN_POINTER(result);
3281 }
3282 
3283 /*
3284  * ST_Snap
3285  *
3286  * Snap a geometry to another with a given tolerance
3287  */
3288 Datum ST_Snap(PG_FUNCTION_ARGS);
3289 PG_FUNCTION_INFO_V1(ST_Snap);
3290 Datum ST_Snap(PG_FUNCTION_ARGS)
3291 {
3292  GSERIALIZED *geom1, *geom2, *result;
3293  LWGEOM *lwgeom1, *lwgeom2, *lwresult;
3294  double tolerance;
3295 
3296  geom1 = PG_GETARG_GSERIALIZED_P(0);
3297  geom2 = PG_GETARG_GSERIALIZED_P(1);
3298  tolerance = PG_GETARG_FLOAT8(2);
3299 
3300  lwgeom1 = lwgeom_from_gserialized(geom1);
3301  lwgeom2 = lwgeom_from_gserialized(geom2);
3302 
3303  lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3304  lwgeom_free(lwgeom1);
3305  lwgeom_free(lwgeom2);
3306  PG_FREE_IF_COPY(geom1, 0);
3307  PG_FREE_IF_COPY(geom2, 1);
3308 
3309  result = geometry_serialize(lwresult);
3310  lwgeom_free(lwresult);
3311 
3312  PG_RETURN_POINTER(result);
3313 }
3314 
3315 /*
3316  * ST_Split
3317  *
3318  * Split polygon by line, line by line, line by point.
3319  * Returns at most components as a collection.
3320  * First element of the collection is always the part which
3321  * remains after the cut, while the second element is the
3322  * part which has been cut out. We arbitrarily take the part
3323  * on the *right* of cut lines as the part which has been cut out.
3324  * For a line cut by a point the part which remains is the one
3325  * from start of the line to the cut point.
3326  *
3327  *
3328  * Author: Sandro Santilli <strk@kbt.io>
3329  *
3330  * Work done for Faunalia (http://www.faunalia.it) with fundings
3331  * from Regione Toscana - Sistema Informativo per il Governo
3332  * del Territorio e dell'Ambiente (RT-SIGTA).
3333  *
3334  * Thanks to the PostGIS community for sharing poly/line ideas [1]
3335  *
3336  * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3337  *
3338  */
3339 Datum ST_Split(PG_FUNCTION_ARGS);
3340 PG_FUNCTION_INFO_V1(ST_Split);
3341 Datum ST_Split(PG_FUNCTION_ARGS)
3342 {
3343  GSERIALIZED *in, *blade_in, *out;
3344  LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
3345 
3346  in = PG_GETARG_GSERIALIZED_P(0);
3347  blade_in = PG_GETARG_GSERIALIZED_P(1);
3348  gserialized_error_if_srid_mismatch(in, blade_in, __func__);
3349 
3350  lwgeom_in = lwgeom_from_gserialized(in);
3351  lwblade_in = lwgeom_from_gserialized(blade_in);
3352 
3353  if (!lwgeom_isfinite(lwgeom_in))
3354  {
3355  lwpgerror("Input Geometry contains invalid coordinates");
3356  PG_RETURN_NULL();
3357  }
3358 
3359  if (!lwgeom_isfinite(lwblade_in))
3360  {
3361  lwpgerror("Blade Geometry contains invalid coordinates");
3362  PG_RETURN_NULL();
3363  }
3364 
3365 
3366  lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3367  lwgeom_free(lwgeom_in);
3368  lwgeom_free(lwblade_in);
3369 
3370  if ( ! lwgeom_out )
3371  {
3372  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3373  PG_FREE_IF_COPY(blade_in, 1);
3374  PG_RETURN_NULL();
3375  }
3376 
3377  out = geometry_serialize(lwgeom_out);
3378  lwgeom_free(lwgeom_out);
3379  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3380  PG_FREE_IF_COPY(blade_in, 1);
3381 
3382  PG_RETURN_POINTER(out);
3383 }
3384 
3385 /**********************************************************************
3386  *
3387  * ST_SharedPaths
3388  *
3389  * Return the set of paths shared between two linear geometries,
3390  * and their direction (same or opposite).
3391  *
3392  * Developed by Sandro Santilli (strk@kbt.io) for Faunalia
3393  * (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3394  * Informativo per la Gestione del Territorio e dell' Ambiente
3395  * [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3396  * trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3397  * 0494241492)"
3398  *
3399  **********************************************************************/
3400 Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3401 PG_FUNCTION_INFO_V1(ST_SharedPaths);
3402 Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3403 {
3404  GSERIALIZED *geom1, *geom2, *out;
3405  LWGEOM *g1, *g2, *lwgeom_out;
3406 
3407  geom1 = PG_GETARG_GSERIALIZED_P(0);
3408  geom2 = PG_GETARG_GSERIALIZED_P(1);
3409 
3410  g1 = lwgeom_from_gserialized(geom1);
3411  g2 = lwgeom_from_gserialized(geom2);
3412 
3413  lwgeom_out = lwgeom_sharedpaths(g1, g2);
3414  lwgeom_free(g1);
3415  lwgeom_free(g2);
3416 
3417  if ( ! lwgeom_out )
3418  {
3419  PG_FREE_IF_COPY(geom1, 0);
3420  PG_FREE_IF_COPY(geom2, 1);
3421  PG_RETURN_NULL();
3422  }
3423 
3424  out = geometry_serialize(lwgeom_out);
3425  lwgeom_free(lwgeom_out);
3426 
3427  PG_FREE_IF_COPY(geom1, 0);
3428  PG_FREE_IF_COPY(geom2, 1);
3429  PG_RETURN_POINTER(out);
3430 }
3431 
3432 
3433 /**********************************************************************
3434  *
3435  * ST_Node
3436  *
3437  * Fully node a set of lines using the least possible nodes while
3438  * preserving all of the input ones.
3439  *
3440  **********************************************************************/
3441 Datum ST_Node(PG_FUNCTION_ARGS);
3442 PG_FUNCTION_INFO_V1(ST_Node);
3443 Datum ST_Node(PG_FUNCTION_ARGS)
3444 {
3445  GSERIALIZED *geom1, *out;
3446  LWGEOM *g1, *lwgeom_out;
3447 
3448  geom1 = PG_GETARG_GSERIALIZED_P(0);
3449 
3450  g1 = lwgeom_from_gserialized(geom1);
3451 
3452  lwgeom_out = lwgeom_node(g1);
3453  lwgeom_free(g1);
3454 
3455  if ( ! lwgeom_out )
3456  {
3457  PG_FREE_IF_COPY(geom1, 0);
3458  PG_RETURN_NULL();
3459  }
3460 
3461  out = geometry_serialize(lwgeom_out);
3462  lwgeom_free(lwgeom_out);
3463 
3464  PG_FREE_IF_COPY(geom1, 0);
3465  PG_RETURN_POINTER(out);
3466 }
3467 
3468 /******************************************
3469  *
3470  * ST_Voronoi
3471  *
3472  * Returns a Voronoi diagram constructed
3473  * from the points of the input geometry.
3474  *
3475  ******************************************/
3476 Datum ST_Voronoi(PG_FUNCTION_ARGS);
3477 PG_FUNCTION_INFO_V1(ST_Voronoi);
3478 Datum ST_Voronoi(PG_FUNCTION_ARGS)
3479 {
3480  GSERIALIZED* input;
3481  GSERIALIZED* clip;
3482  GSERIALIZED* result;
3483  LWGEOM* lwgeom_input;
3484  LWGEOM* lwgeom_result;
3485  double tolerance;
3486  GBOX clip_envelope;
3487  int custom_clip_envelope;
3488  int return_polygons;
3489 
3490  /* Return NULL on NULL geometry */
3491  if (PG_ARGISNULL(0))
3492  PG_RETURN_NULL();
3493 
3494  /* Read our tolerance value */
3495  if (PG_ARGISNULL(2))
3496  {
3497  lwpgerror("Tolerance must be a positive number.");
3498  PG_RETURN_NULL();
3499  }
3500 
3501  tolerance = PG_GETARG_FLOAT8(2);
3502 
3503  if (tolerance < 0)
3504  {
3505  lwpgerror("Tolerance must be a positive number.");
3506  PG_RETURN_NULL();
3507  }
3508 
3509  /* Are we returning lines or polygons? */
3510  if (PG_ARGISNULL(3))
3511  {
3512  lwpgerror("return_polygons must be true or false.");
3513  PG_RETURN_NULL();
3514  }
3515  return_polygons = PG_GETARG_BOOL(3);
3516 
3517  /* Read our clipping envelope, if applicable. */
3518  custom_clip_envelope = !PG_ARGISNULL(1);
3519  if (custom_clip_envelope) {
3520  clip = PG_GETARG_GSERIALIZED_P(1);
3521  if (!gserialized_get_gbox_p(clip, &clip_envelope))
3522  {
3523  lwpgerror("Could not determine envelope of clipping geometry.");
3524  PG_FREE_IF_COPY(clip, 1);
3525  PG_RETURN_NULL();
3526  }
3527  PG_FREE_IF_COPY(clip, 1);
3528  }
3529 
3530  /* Read our input geometry */
3531  input = PG_GETARG_GSERIALIZED_P(0);
3532 
3533  lwgeom_input = lwgeom_from_gserialized(input);
3534 
3535  if(!lwgeom_input)
3536  {
3537  lwpgerror("Could not read input geometry.");
3538  PG_FREE_IF_COPY(input, 0);
3539  PG_RETURN_NULL();
3540  }
3541 
3542  lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
3543  lwgeom_free(lwgeom_input);
3544 
3545  if (!lwgeom_result)
3546  {
3547  lwpgerror("Error computing Voronoi diagram.");
3548  PG_FREE_IF_COPY(input, 0);
3549  PG_RETURN_NULL();
3550  }
3551 
3552  result = geometry_serialize(lwgeom_result);
3553  lwgeom_free(lwgeom_result);
3554 
3555  PG_FREE_IF_COPY(input, 0);
3556  PG_RETURN_POINTER(result);
3557 }
3558 
3559 /******************************************
3560  *
3561  * ST_MinimumClearance
3562  *
3563  * Returns the minimum clearance of a geometry.
3564  *
3565  ******************************************/
3566 Datum ST_MinimumClearance(PG_FUNCTION_ARGS);
3567 PG_FUNCTION_INFO_V1(ST_MinimumClearance);
3568 Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
3569 {
3570  GSERIALIZED* input;
3571  GEOSGeometry* input_geos;
3572  int error;
3573  double result;
3574 
3575  initGEOS(lwpgnotice, lwgeom_geos_error);
3576 
3577  input = PG_GETARG_GSERIALIZED_P(0);
3578  input_geos = POSTGIS2GEOS(input);
3579  if (!input_geos)
3580  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3581 
3582  error = GEOSMinimumClearance(input_geos, &result);
3583  GEOSGeom_destroy(input_geos);
3584  if (error) HANDLE_GEOS_ERROR("Error computing minimum clearance");
3585 
3586  PG_FREE_IF_COPY(input, 0);
3587  PG_RETURN_FLOAT8(result);
3588 }
3589 
3590 /******************************************
3591  *
3592  * ST_MinimumClearanceLine
3593  *
3594  * Returns the minimum clearance line of a geometry.
3595  *
3596  ******************************************/
3597 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS);
3598 PG_FUNCTION_INFO_V1(ST_MinimumClearanceLine);
3599 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
3600 {
3601  GSERIALIZED* input;
3602  GSERIALIZED* result;
3603  GEOSGeometry* input_geos;
3604  GEOSGeometry* result_geos;
3605  int32_t srid;
3606 
3607  initGEOS(lwpgnotice, lwgeom_geos_error);
3608 
3609  input = PG_GETARG_GSERIALIZED_P(0);
3610  srid = gserialized_get_srid(input);
3611  input_geos = POSTGIS2GEOS(input);
3612  if (!input_geos)
3613  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3614 
3615  result_geos = GEOSMinimumClearanceLine(input_geos);
3616  GEOSGeom_destroy(input_geos);
3617  if (!result_geos)
3618  HANDLE_GEOS_ERROR("Error computing minimum clearance");
3619 
3620  GEOSSetSRID(result_geos, srid);
3621  result = GEOS2POSTGIS(result_geos, LW_FALSE);
3622  GEOSGeom_destroy(result_geos);
3623 
3624  PG_FREE_IF_COPY(input, 0);
3625  PG_RETURN_POINTER(result);
3626 }
3627 
3628 /******************************************
3629  *
3630  * ST_OrientedEnvelope
3631  *
3632  ******************************************/
3633 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS);
3634 PG_FUNCTION_INFO_V1(ST_OrientedEnvelope);
3635 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
3636 {
3637  GSERIALIZED* input;
3638  GSERIALIZED* result;
3639  GEOSGeometry* input_geos;
3640  GEOSGeometry* result_geos;
3641  int32_t srid;
3642 
3643  initGEOS(lwpgnotice, lwgeom_geos_error);
3644 
3645  input = PG_GETARG_GSERIALIZED_P(0);
3646  srid = gserialized_get_srid(input);
3647  input_geos = POSTGIS2GEOS(input);
3648  if (!input_geos)
3649  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3650 
3651  result_geos = GEOSMinimumRotatedRectangle(input_geos);
3652  GEOSGeom_destroy(input_geos);
3653  if (!result_geos)
3654  HANDLE_GEOS_ERROR("Error computing oriented envelope");
3655 
3656  GEOSSetSRID(result_geos, srid);
3657  result = GEOS2POSTGIS(result_geos, LW_FALSE);
3658  GEOSGeom_destroy(result_geos);
3659 
3660  PG_FREE_IF_COPY(input, 0);
3661  PG_RETURN_POINTER(result);
3662 }
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:2859
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:2994
contains
Datum contains(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:1842
coveredby
Datum coveredby(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2180
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:1696
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:1620
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:3290
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:3224
ST_Split
Datum ST_Split(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3341
covers
Datum covers(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2049
ST_Node
Datum ST_Node(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3443
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:3118
centroid
Datum centroid(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:1522
touches
Datum touches(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2469
disjoint
Datum disjoint(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2526
ST_MinimumClearance
Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3568
ST_Intersects
Datum ST_Intersects(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2348
GEOSnoop
Datum GEOSnoop(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2973
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:2292
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:2641
relate_pattern
Datum relate_pattern(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2583
issimple
Datum issimple(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2759
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:2875
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:2819
ST_DelaunayTriangles
Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3254
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:2844
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:1785
ST_Equals
Datum ST_Equals(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2693
within
Datum within(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:1974
ST_OrientedEnvelope
Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3635
ST_SharedPaths
Datum ST_SharedPaths(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3402
ST_MinimumClearanceLine
Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3599
ST_Voronoi
Datum ST_Voronoi(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3478
clusterintersecting_garray
Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3052
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:2917
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:1666
isring
Datum isring(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:2785
linemerge
Datum linemerge(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:3194
pointonsurface
Datum pointonsurface(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:1502
containsproperly
Datum containsproperly(PG_FUNCTION_ARGS)
Definition: postgis/lwgeom_geos.c:1983
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