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