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