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