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