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