PostGIS  2.5.0beta2dev-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 
2152  errorIfGeometryCollection(geom1,geom2);
2154 
2155  /* A.Intersects(Empty) == FALSE */
2156  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2157  PG_RETURN_BOOL(false);
2158 
2159  /*
2160  * short-circuit 1: if geom2 bounding box does not overlap
2161  * geom1 bounding box we can return FALSE.
2162  */
2163  if ( gserialized_get_gbox_p(geom1, &box1) &&
2164  gserialized_get_gbox_p(geom2, &box2) )
2165  {
2166  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2167  {
2168  PG_RETURN_BOOL(false);
2169  }
2170  }
2171 
2172  /*
2173  * short-circuit 2: if the geoms are a point and a polygon,
2174  * call the point_outside_polygon function.
2175  */
2176  if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
2177  {
2178  GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
2179  GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
2180  RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
2181  int retval;
2182 
2183  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2184  if (gserialized_get_type(gpoint) == POINTTYPE)
2185  {
2186  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2187  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2188  lwgeom_free(point);
2189 
2190  retval = (pip_result != -1); /* not outside */
2191  }
2192  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2193  {
2195  uint32_t i;
2196 
2197  retval = LW_FALSE;
2198  for (i = 0; i < mpoint->ngeoms; i++)
2199  {
2200  int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2201  if (pip_result != -1) /* not outside */
2202  {
2203  retval = LW_TRUE;
2204  break;
2205  }
2206  }
2207 
2208  lwmpoint_free(mpoint);
2209  }
2210  else
2211  {
2212  /* Never get here */
2213  elog(ERROR,"Type isn't point or multipoint!");
2214  PG_RETURN_NULL();
2215  }
2216 
2217  PG_FREE_IF_COPY(geom1, 0);
2218  PG_FREE_IF_COPY(geom2, 1);
2219  PG_RETURN_BOOL(retval);
2220  }
2221 
2222  initGEOS(lwpgnotice, lwgeom_geos_error);
2223  prep_cache = GetPrepGeomCache( fcinfo, geom1, geom2 );
2224 
2225  if ( prep_cache && prep_cache->prepared_geom )
2226  {
2227  if ( prep_cache->gcache.argnum == 1 )
2228  {
2229  GEOSGeometry *g = POSTGIS2GEOS(geom2);
2230  if (!g) HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2231  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2232  GEOSGeom_destroy(g);
2233  }
2234  else
2235  {
2236  GEOSGeometry *g = POSTGIS2GEOS(geom1);
2237  if (!g)
2238  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2239  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2240  GEOSGeom_destroy(g);
2241  }
2242  }
2243  else
2244  {
2245  GEOSGeometry *g1;
2246  GEOSGeometry *g2;
2247  g1 = POSTGIS2GEOS(geom1);
2248  if (!g1) HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2249  g2 = POSTGIS2GEOS(geom2);
2250  if (!g2)
2251  {
2252  GEOSGeom_destroy(g1);
2253  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2254  }
2255  result = GEOSIntersects( g1, g2);
2256  GEOSGeom_destroy(g1);
2257  GEOSGeom_destroy(g2);
2258  }
2259 
2260  if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
2261 
2262  PG_FREE_IF_COPY(geom1, 0);
2263  PG_FREE_IF_COPY(geom2, 1);
2264 
2265  PG_RETURN_BOOL(result);
2266 }
2267 
2268 
2270 Datum touches(PG_FUNCTION_ARGS)
2271 {
2272  GSERIALIZED *geom1;
2273  GSERIALIZED *geom2;
2274  GEOSGeometry *g1, *g2;
2275  char result;
2276  GBOX box1, box2;
2277 
2278  geom1 = PG_GETARG_GSERIALIZED_P(0);
2279  geom2 = PG_GETARG_GSERIALIZED_P(1);
2280 
2281  errorIfGeometryCollection(geom1,geom2);
2283 
2284  /* A.Touches(Empty) == FALSE */
2285  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2286  PG_RETURN_BOOL(false);
2287 
2288  /*
2289  * short-circuit 1: if geom2 bounding box does not overlap
2290  * geom1 bounding box we can return FALSE.
2291  */
2292  if ( gserialized_get_gbox_p(geom1, &box1) &&
2293  gserialized_get_gbox_p(geom2, &box2) )
2294  {
2295  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2296  {
2297  PG_RETURN_BOOL(false);
2298  }
2299  }
2300 
2301  initGEOS(lwpgnotice, lwgeom_geos_error);
2302 
2303  g1 = POSTGIS2GEOS(geom1 );
2304  if (!g1)
2305  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2306 
2307  g2 = POSTGIS2GEOS(geom2 );
2308  if (!g2)
2309  {
2310  GEOSGeom_destroy(g1);
2311  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2312  }
2313 
2314  result = GEOSTouches(g1,g2);
2315 
2316  GEOSGeom_destroy(g1);
2317  GEOSGeom_destroy(g2);
2318 
2319  if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
2320 
2321  PG_FREE_IF_COPY(geom1, 0);
2322  PG_FREE_IF_COPY(geom2, 1);
2323 
2324  PG_RETURN_BOOL(result);
2325 }
2326 
2327 
2329 Datum disjoint(PG_FUNCTION_ARGS)
2330 {
2331  GSERIALIZED *geom1;
2332  GSERIALIZED *geom2;
2333  GEOSGeometry *g1, *g2;
2334  char result;
2335  GBOX box1, box2;
2336 
2337  geom1 = PG_GETARG_GSERIALIZED_P(0);
2338  geom2 = PG_GETARG_GSERIALIZED_P(1);
2339 
2340  errorIfGeometryCollection(geom1,geom2);
2342 
2343  /* A.Disjoint(Empty) == TRUE */
2344  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2345  PG_RETURN_BOOL(true);
2346 
2347  /*
2348  * short-circuit 1: if geom2 bounding box does not overlap
2349  * geom1 bounding box we can return TRUE.
2350  */
2351  if ( gserialized_get_gbox_p(geom1, &box1) &&
2352  gserialized_get_gbox_p(geom2, &box2) )
2353  {
2354  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2355  {
2356  PG_RETURN_BOOL(true);
2357  }
2358  }
2359 
2360  initGEOS(lwpgnotice, lwgeom_geos_error);
2361 
2362  g1 = POSTGIS2GEOS(geom1);
2363  if (!g1)
2364  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2365 
2366  g2 = POSTGIS2GEOS(geom2);
2367  if (!g2)
2368  {
2369  GEOSGeom_destroy(g1);
2370  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2371  }
2372 
2373  result = GEOSDisjoint(g1,g2);
2374 
2375  GEOSGeom_destroy(g1);
2376  GEOSGeom_destroy(g2);
2377 
2378  if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
2379 
2380  PG_FREE_IF_COPY(geom1, 0);
2381  PG_FREE_IF_COPY(geom2, 1);
2382 
2383  PG_RETURN_BOOL(result);
2384 }
2385 
2386 
2388 Datum relate_pattern(PG_FUNCTION_ARGS)
2389 {
2390  GSERIALIZED *geom1;
2391  GSERIALIZED *geom2;
2392  char *patt;
2393  char result;
2394  GEOSGeometry *g1, *g2;
2395  size_t i;
2396 
2397  geom1 = PG_GETARG_GSERIALIZED_P(0);
2398  geom2 = PG_GETARG_GSERIALIZED_P(1);
2399 
2400  /* TODO handle empty */
2401 
2402  errorIfGeometryCollection(geom1,geom2);
2404 
2405  initGEOS(lwpgnotice, lwgeom_geos_error);
2406 
2407  g1 = POSTGIS2GEOS(geom1);
2408  if (!g1)
2409  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2410  g2 = POSTGIS2GEOS(geom2);
2411  if (!g2)
2412  {
2413  GEOSGeom_destroy(g1);
2414  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2415  }
2416 
2417  patt = DatumGetCString(DirectFunctionCall1(textout,
2418  PointerGetDatum(PG_GETARG_DATUM(2))));
2419 
2420  /*
2421  ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2422  */
2423  for ( i = 0; i < strlen(patt); i++ )
2424  {
2425  if ( patt[i] == 't' ) patt[i] = 'T';
2426  if ( patt[i] == 'f' ) patt[i] = 'F';
2427  }
2428 
2429  result = GEOSRelatePattern(g1,g2,patt);
2430  GEOSGeom_destroy(g1);
2431  GEOSGeom_destroy(g2);
2432  pfree(patt);
2433 
2434  if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
2435 
2436  PG_FREE_IF_COPY(geom1, 0);
2437  PG_FREE_IF_COPY(geom2, 1);
2438 
2439  PG_RETURN_BOOL(result);
2440 }
2441 
2442 
2443 
2445 Datum relate_full(PG_FUNCTION_ARGS)
2446 {
2447  GSERIALIZED *geom1;
2448  GSERIALIZED *geom2;
2449  GEOSGeometry *g1, *g2;
2450  char *relate_str;
2451  text *result;
2452  int bnr = GEOSRELATE_BNR_OGC;
2453 
2454  POSTGIS_DEBUG(2, "in relate_full()");
2455 
2456  /* TODO handle empty */
2457 
2458  geom1 = PG_GETARG_GSERIALIZED_P(0);
2459  geom2 = PG_GETARG_GSERIALIZED_P(1);
2460 
2461  if ( PG_NARGS() > 2 )
2462  {
2463  bnr = PG_GETARG_INT32(2);
2464  }
2465 
2466  errorIfGeometryCollection(geom1,geom2);
2468 
2469  initGEOS(lwpgnotice, lwgeom_geos_error);
2470 
2471  g1 = POSTGIS2GEOS(geom1 );
2472  if (!g1)
2473  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2474  g2 = POSTGIS2GEOS(geom2 );
2475  if (!g2)
2476  {
2477  GEOSGeom_destroy(g1);
2478  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2479  }
2480 
2481  POSTGIS_DEBUG(3, "constructed geometries ");
2482 
2483  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2484  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2485 
2486  relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2487 
2488  GEOSGeom_destroy(g1);
2489  GEOSGeom_destroy(g2);
2490 
2491  if (!relate_str) HANDLE_GEOS_ERROR("GEOSRelate");
2492 
2493  result = cstring_to_text(relate_str);
2494  GEOSFree(relate_str);
2495 
2496  PG_FREE_IF_COPY(geom1, 0);
2497  PG_FREE_IF_COPY(geom2, 1);
2498 
2499  PG_RETURN_TEXT_P(result);
2500 }
2501 
2502 
2504 Datum ST_Equals(PG_FUNCTION_ARGS)
2505 {
2506  GSERIALIZED *geom1;
2507  GSERIALIZED *geom2;
2508  GEOSGeometry *g1, *g2;
2509  char result;
2510  GBOX box1, box2;
2511 
2512  geom1 = PG_GETARG_GSERIALIZED_P(0);
2513  geom2 = PG_GETARG_GSERIALIZED_P(1);
2514 
2515  errorIfGeometryCollection(geom1,geom2);
2517 
2518  /* Empty == Empty */
2519  if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
2520  PG_RETURN_BOOL(true);
2521 
2522  /*
2523  * short-circuit: If geom1 and geom2 do not have the same bounding box
2524  * we can return FALSE.
2525  */
2526  if ( gserialized_get_gbox_p(geom1, &box1) &&
2527  gserialized_get_gbox_p(geom2, &box2) )
2528  {
2529  if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
2530  {
2531  PG_RETURN_BOOL(false);
2532  }
2533  }
2534 
2535  /*
2536  * short-circuit: if geom1 and geom2 are binary-equivalent, we can return
2537  * TRUE. This is much faster than doing the comparison using GEOS.
2538  */
2539  if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
2540  PG_RETURN_BOOL(true);
2541  }
2542 
2543  initGEOS(lwpgnotice, lwgeom_geos_error);
2544 
2545  g1 = POSTGIS2GEOS(geom1);
2546 
2547  if (!g1)
2548  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2549 
2550  g2 = POSTGIS2GEOS(geom2);
2551 
2552  if (!g2)
2553  {
2554  GEOSGeom_destroy(g1);
2555  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
2556  }
2557 
2558  result = GEOSEquals(g1,g2);
2559 
2560  GEOSGeom_destroy(g1);
2561  GEOSGeom_destroy(g2);
2562 
2563  if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
2564 
2565  PG_FREE_IF_COPY(geom1, 0);
2566  PG_FREE_IF_COPY(geom2, 1);
2567 
2568  PG_RETURN_BOOL(result);
2569 }
2570 
2572 Datum issimple(PG_FUNCTION_ARGS)
2573 {
2574  GSERIALIZED *geom;
2575  LWGEOM *lwgeom_in;
2576  int result;
2577 
2578  POSTGIS_DEBUG(2, "issimple called");
2579 
2580  geom = PG_GETARG_GSERIALIZED_P(0);
2581 
2582  if ( gserialized_is_empty(geom) )
2583  PG_RETURN_BOOL(true);
2584 
2585  lwgeom_in = lwgeom_from_gserialized(geom);
2586  result = lwgeom_is_simple(lwgeom_in);
2587  lwgeom_free(lwgeom_in) ;
2588  PG_FREE_IF_COPY(geom, 0);
2589 
2590  if (result == -1) {
2591  PG_RETURN_NULL(); /*never get here */
2592  }
2593 
2594  PG_RETURN_BOOL(result);
2595 }
2596 
2598 Datum isring(PG_FUNCTION_ARGS)
2599 {
2600  GSERIALIZED *geom;
2601  GEOSGeometry *g1;
2602  int result;
2603 
2604  geom = PG_GETARG_GSERIALIZED_P(0);
2605 
2606  /* Empty things can't close */
2607  if ( gserialized_is_empty(geom) )
2608  PG_RETURN_BOOL(false);
2609 
2610  initGEOS(lwpgnotice, lwgeom_geos_error);
2611 
2612  g1 = POSTGIS2GEOS(geom);
2613  if (!g1)
2614  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2615 
2616  if ( GEOSGeomTypeId(g1) != GEOS_LINESTRING )
2617  {
2618  GEOSGeom_destroy(g1);
2619  elog(ERROR, "ST_IsRing() should only be called on a linear feature");
2620  }
2621 
2622  result = GEOSisRing(g1);
2623  GEOSGeom_destroy(g1);
2624 
2625  if (result == 2) HANDLE_GEOS_ERROR("GEOSisRing");
2626 
2627  PG_FREE_IF_COPY(geom, 0);
2628  PG_RETURN_BOOL(result);
2629 }
2630 
2631 GSERIALIZED *
2632 GEOS2POSTGIS(GEOSGeom geom, char want3d)
2633 {
2634  LWGEOM *lwgeom;
2635  GSERIALIZED *result;
2636 
2637  lwgeom = GEOS2LWGEOM(geom, want3d);
2638  if ( ! lwgeom )
2639  {
2640  lwpgerror("%s: GEOS2LWGEOM returned NULL", __func__);
2641  return NULL;
2642  }
2643 
2644  POSTGIS_DEBUGF(4, "%s: GEOS2LWGEOM returned a %s", __func__, lwgeom_summary(lwgeom, 0));
2645 
2646  if (lwgeom_needs_bbox(lwgeom)) lwgeom_add_bbox(lwgeom);
2647 
2648  result = geometry_serialize(lwgeom);
2649  lwgeom_free(lwgeom);
2650 
2651  return result;
2652 }
2653 
2654 /*-----=POSTGIS2GEOS= */
2655 
2656 
2657 GEOSGeometry *
2659 {
2660  GEOSGeometry *ret;
2661  LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
2662  if ( ! lwgeom )
2663  {
2664  lwpgerror("POSTGIS2GEOS: unable to deserialize input");
2665  return NULL;
2666  }
2667  ret = LWGEOM2GEOS(lwgeom, 0);
2668  lwgeom_free(lwgeom);
2669 
2670  return ret;
2671 }
2672 
2673 uint32_t array_nelems_not_null(ArrayType* array) {
2674  ArrayIterator iterator;
2675  Datum value;
2676  bool isnull;
2677  uint32_t nelems_not_null = 0;
2678 
2679 #if POSTGIS_PGSQL_VERSION >= 95
2680  iterator = array_create_iterator(array, 0, NULL);
2681 #else
2682  iterator = array_create_iterator(array, 0);
2683 #endif
2684  while(array_iterate(iterator, &value, &isnull) )
2685  {
2686  if (!isnull)
2687  {
2688  nelems_not_null++;
2689  }
2690  }
2691  array_free_iterator(iterator);
2692 
2693  return nelems_not_null;
2694 }
2695 
2696 /* ARRAY2LWGEOM: Converts the non-null elements of a Postgres array into a LWGEOM* array */
2697 LWGEOM** ARRAY2LWGEOM(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2698 {
2699  ArrayIterator iterator;
2700  Datum value;
2701  bool isnull;
2702  bool gotsrid = false;
2703  uint32_t i = 0;
2704 
2705  LWGEOM** lw_geoms = palloc(nelems * sizeof(LWGEOM*));
2706 
2707 #if POSTGIS_PGSQL_VERSION >= 95
2708  iterator = array_create_iterator(array, 0, NULL);
2709 #else
2710  iterator = array_create_iterator(array, 0);
2711 #endif
2712 
2713  while(array_iterate(iterator, &value, &isnull))
2714  {
2715  GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2716 
2717  if (isnull)
2718  {
2719  continue;
2720  }
2721 
2722  *is3d = *is3d || gserialized_has_z(geom);
2723 
2724  lw_geoms[i] = lwgeom_from_gserialized(geom);
2725  if (!lw_geoms[i]) /* error in creation */
2726  {
2727  lwpgerror("Geometry deserializing geometry");
2728  return NULL;
2729  }
2730  if (!gotsrid)
2731  {
2732  gotsrid = true;
2733  *srid = gserialized_get_srid(geom);
2734  }
2735  else if (*srid != gserialized_get_srid(geom))
2736  {
2738  return NULL;
2739  }
2740 
2741  i++;
2742  }
2743 
2744  return lw_geoms;
2745 }
2746 
2747 /* ARRAY2GEOS: Converts the non-null elements of a Postgres array into a GEOSGeometry* array */
2748 GEOSGeometry** ARRAY2GEOS(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2749 {
2750  ArrayIterator iterator;
2751  Datum value;
2752  bool isnull;
2753  bool gotsrid = false;
2754  uint32_t i = 0;
2755 
2756  GEOSGeometry** geos_geoms = palloc(nelems * sizeof(GEOSGeometry*));
2757 
2758 #if POSTGIS_PGSQL_VERSION >= 95
2759  iterator = array_create_iterator(array, 0, NULL);
2760 #else
2761  iterator = array_create_iterator(array, 0);
2762 #endif
2763 
2764  while(array_iterate(iterator, &value, &isnull))
2765  {
2766  GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2767 
2768  if (isnull)
2769  {
2770  continue;
2771  }
2772 
2773  *is3d = *is3d || gserialized_has_z(geom);
2774 
2775  geos_geoms[i] = POSTGIS2GEOS(geom);
2776  if (!geos_geoms[i])
2777  {
2778  uint32_t j;
2779  lwpgerror("Geometry could not be converted to GEOS");
2780 
2781  for (j = 0; j < i; j++) {
2782  GEOSGeom_destroy(geos_geoms[j]);
2783  }
2784  return NULL;
2785  }
2786 
2787  if (!gotsrid)
2788  {
2789  *srid = gserialized_get_srid(geom);
2790  gotsrid = true;
2791  }
2792  else if (*srid != gserialized_get_srid(geom))
2793  {
2794  uint32_t j;
2796 
2797  for (j = 0; j <= i; j++) {
2798  GEOSGeom_destroy(geos_geoms[j]);
2799  }
2800  return NULL;
2801  }
2802 
2803  i++;
2804  }
2805 
2806  array_free_iterator(iterator);
2807  return geos_geoms;
2808 }
2809 
2811 Datum GEOSnoop(PG_FUNCTION_ARGS)
2812 {
2813  GSERIALIZED *geom;
2814  GEOSGeometry *geosgeom;
2815  GSERIALIZED *lwgeom_result;
2816 
2817  initGEOS(lwpgnotice, lwgeom_geos_error);
2818 
2819  geom = PG_GETARG_GSERIALIZED_P(0);
2820  geosgeom = POSTGIS2GEOS(geom);
2821  if ( ! geosgeom ) PG_RETURN_NULL();
2822 
2823  lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
2824  GEOSGeom_destroy(geosgeom);
2825 
2826  PG_FREE_IF_COPY(geom, 0);
2827 
2828  PG_RETURN_POINTER(lwgeom_result);
2829 }
2830 
2832 Datum polygonize_garray(PG_FUNCTION_ARGS)
2833 {
2834  ArrayType *array;
2835  int is3d = 0;
2836  uint32 nelems, i;
2837  GSERIALIZED *result;
2838  GEOSGeometry *geos_result;
2839  const GEOSGeometry **vgeoms;
2840  int srid=SRID_UNKNOWN;
2841 #if POSTGIS_DEBUG_LEVEL >= 3
2842  static int call=1;
2843 #endif
2844 
2845 #if POSTGIS_DEBUG_LEVEL >= 3
2846  call++;
2847 #endif
2848 
2849  if (PG_ARGISNULL(0))
2850  PG_RETURN_NULL();
2851 
2852  array = PG_GETARG_ARRAYTYPE_P(0);
2853  nelems = array_nelems_not_null(array);
2854 
2855  if (nelems == 0)
2856  PG_RETURN_NULL();
2857 
2858  POSTGIS_DEBUGF(3, "polygonize_garray: number of non-null elements: %d", nelems);
2859 
2860  /* Ok, we really need geos now ;) */
2861  initGEOS(lwpgnotice, lwgeom_geos_error);
2862 
2863  vgeoms = (const GEOSGeometry**) ARRAY2GEOS(array, nelems, &is3d, &srid);
2864 
2865  POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
2866 
2867  geos_result = GEOSPolygonize(vgeoms, nelems);
2868 
2869  POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
2870 
2871  for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
2872  pfree(vgeoms);
2873 
2874  if ( ! geos_result ) PG_RETURN_NULL();
2875 
2876  GEOSSetSRID(geos_result, srid);
2877  result = GEOS2POSTGIS(geos_result, is3d);
2878  GEOSGeom_destroy(geos_result);
2879  if (!result)
2880  {
2881  elog(ERROR, "%s returned an error", __func__);
2882  PG_RETURN_NULL(); /*never get here */
2883  }
2884 
2885  PG_RETURN_POINTER(result);
2886 }
2887 
2888 
2890 Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
2891 {
2892  Datum* result_array_data;
2893  ArrayType *array, *result;
2894  int is3d = 0;
2895  uint32 nelems, nclusters, i;
2896  GEOSGeometry **geos_inputs, **geos_results;
2897  int srid=SRID_UNKNOWN;
2898 
2899  /* Parameters used to construct a result array */
2900  int16 elmlen;
2901  bool elmbyval;
2902  char elmalign;
2903 
2904  /* Null array, null geometry (should be empty?) */
2905  if (PG_ARGISNULL(0))
2906  PG_RETURN_NULL();
2907 
2908  array = PG_GETARG_ARRAYTYPE_P(0);
2909  nelems = array_nelems_not_null(array);
2910 
2911  POSTGIS_DEBUGF(3, "clusterintersecting_garray: number of non-null elements: %d", nelems);
2912 
2913  if ( nelems == 0 ) PG_RETURN_NULL();
2914 
2915  /* TODO short-circuit for one element? */
2916 
2917  /* Ok, we really need geos now ;) */
2918  initGEOS(lwpgnotice, lwgeom_geos_error);
2919 
2920  geos_inputs = ARRAY2GEOS(array, nelems, &is3d, &srid);
2921  if(!geos_inputs)
2922  {
2923  PG_RETURN_NULL();
2924  }
2925 
2926  if (cluster_intersecting(geos_inputs, nelems, &geos_results, &nclusters) != LW_SUCCESS)
2927  {
2928  elog(ERROR, "clusterintersecting: Error performing clustering");
2929  PG_RETURN_NULL();
2930  }
2931  pfree(geos_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
2932 
2933  if (!geos_results) PG_RETURN_NULL();
2934 
2935  result_array_data = palloc(nclusters * sizeof(Datum));
2936  for (i=0; i<nclusters; ++i)
2937  {
2938  result_array_data[i] = PointerGetDatum(GEOS2POSTGIS(geos_results[i], is3d));
2939  GEOSGeom_destroy(geos_results[i]);
2940  }
2941  pfree(geos_results);
2942 
2943  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
2944  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
2945 
2946  if (!result)
2947  {
2948  elog(ERROR, "clusterintersecting: Error constructing return-array");
2949  PG_RETURN_NULL();
2950  }
2951 
2952  PG_RETURN_POINTER(result);
2953 }
2954 
2956 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS)
2957 {
2958  Datum* result_array_data;
2959  ArrayType *array, *result;
2960  int is3d = 0;
2961  uint32 nelems, nclusters, i;
2962  LWGEOM** lw_inputs;
2963  LWGEOM** lw_results;
2964  double tolerance;
2965  int srid=SRID_UNKNOWN;
2966 
2967  /* Parameters used to construct a result array */
2968  int16 elmlen;
2969  bool elmbyval;
2970  char elmalign;
2971 
2972  /* Null array, null geometry (should be empty?) */
2973  if (PG_ARGISNULL(0))
2974  PG_RETURN_NULL();
2975 
2976  array = PG_GETARG_ARRAYTYPE_P(0);
2977 
2978  tolerance = PG_GETARG_FLOAT8(1);
2979  if (tolerance < 0)
2980  {
2981  lwpgerror("Tolerance must be a positive number.");
2982  PG_RETURN_NULL();
2983  }
2984 
2985  nelems = array_nelems_not_null(array);
2986 
2987  POSTGIS_DEBUGF(3, "cluster_within_distance_garray: number of non-null elements: %d", nelems);
2988 
2989  if ( nelems == 0 ) PG_RETURN_NULL();
2990 
2991  /* TODO short-circuit for one element? */
2992 
2993  /* Ok, we really need geos now ;) */
2994  initGEOS(lwpgnotice, lwgeom_geos_error);
2995 
2996  lw_inputs = ARRAY2LWGEOM(array, nelems, &is3d, &srid);
2997  if (!lw_inputs)
2998  {
2999  PG_RETURN_NULL();
3000  }
3001 
3002  if (cluster_within_distance(lw_inputs, nelems, tolerance, &lw_results, &nclusters) != LW_SUCCESS)
3003  {
3004  elog(ERROR, "cluster_within: Error performing clustering");
3005  PG_RETURN_NULL();
3006  }
3007  pfree(lw_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3008 
3009  if (!lw_results) PG_RETURN_NULL();
3010 
3011  result_array_data = palloc(nclusters * sizeof(Datum));
3012  for (i=0; i<nclusters; ++i)
3013  {
3014  result_array_data[i] = PointerGetDatum(gserialized_from_lwgeom(lw_results[i], NULL));
3015  lwgeom_free(lw_results[i]);
3016  }
3017  pfree(lw_results);
3018 
3019  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3020  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3021 
3022  if (!result)
3023  {
3024  elog(ERROR, "clusterwithin: Error constructing return-array");
3025  PG_RETURN_NULL();
3026  }
3027 
3028  PG_RETURN_POINTER(result);
3029 }
3030 
3032 Datum linemerge(PG_FUNCTION_ARGS)
3033 {
3034  GSERIALIZED *geom1;
3035  GSERIALIZED *result;
3036  LWGEOM *lwgeom1, *lwresult ;
3037 
3038  geom1 = PG_GETARG_GSERIALIZED_P(0);
3039 
3040 
3041  lwgeom1 = lwgeom_from_gserialized(geom1) ;
3042 
3043  lwresult = lwgeom_linemerge(lwgeom1);
3044  result = geometry_serialize(lwresult) ;
3045 
3046  lwgeom_free(lwgeom1) ;
3047  lwgeom_free(lwresult) ;
3048 
3049  PG_FREE_IF_COPY(geom1, 0);
3050 
3051  PG_RETURN_POINTER(result);
3052 }
3053 
3054 /*
3055  * Take a geometry and return an areal geometry
3056  * (Polygon or MultiPolygon).
3057  * Actually a wrapper around GEOSpolygonize,
3058  * transforming the resulting collection into
3059  * a valid polygon Geometry.
3060  */
3062 Datum ST_BuildArea(PG_FUNCTION_ARGS)
3063 {
3064  GSERIALIZED *result;
3065  GSERIALIZED *geom;
3066  LWGEOM *lwgeom_in, *lwgeom_out;
3067 
3068  geom = PG_GETARG_GSERIALIZED_P(0);
3069  lwgeom_in = lwgeom_from_gserialized(geom);
3070 
3071  lwgeom_out = lwgeom_buildarea(lwgeom_in);
3072  lwgeom_free(lwgeom_in) ;
3073 
3074  if ( ! lwgeom_out )
3075  {
3076  PG_FREE_IF_COPY(geom, 0);
3077  PG_RETURN_NULL();
3078  }
3079 
3080  result = geometry_serialize(lwgeom_out) ;
3081  lwgeom_free(lwgeom_out) ;
3082 
3083  PG_FREE_IF_COPY(geom, 0);
3084  PG_RETURN_POINTER(result);
3085 }
3086 
3087 /*
3088  * Take the vertices of a geometry and builds
3089  * Delaunay triangles around them.
3090  */
3092 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
3093 {
3094  GSERIALIZED *result;
3095  GSERIALIZED *geom;
3096  LWGEOM *lwgeom_in, *lwgeom_out;
3097  double tolerance = 0.0;
3098  int flags = 0;
3099 
3100  geom = PG_GETARG_GSERIALIZED_P(0);
3101  tolerance = PG_GETARG_FLOAT8(1);
3102  flags = PG_GETARG_INT32(2);
3103 
3104  lwgeom_in = lwgeom_from_gserialized(geom);
3105  lwgeom_out = lwgeom_delaunay_triangulation(lwgeom_in, tolerance, flags);
3106  lwgeom_free(lwgeom_in) ;
3107 
3108  if ( ! lwgeom_out )
3109  {
3110  PG_FREE_IF_COPY(geom, 0);
3111  PG_RETURN_NULL();
3112  }
3113 
3114  result = geometry_serialize(lwgeom_out) ;
3115  lwgeom_free(lwgeom_out) ;
3116 
3117  PG_FREE_IF_COPY(geom, 0);
3118  PG_RETURN_POINTER(result);
3119 }
3120 
3121 /*
3122  * ST_Snap
3123  *
3124  * Snap a geometry to another with a given tolerance
3125  */
3126 Datum ST_Snap(PG_FUNCTION_ARGS);
3128 Datum ST_Snap(PG_FUNCTION_ARGS)
3129 {
3130  GSERIALIZED *geom1, *geom2, *result;
3131  LWGEOM *lwgeom1, *lwgeom2, *lwresult;
3132  double tolerance;
3133 
3134  geom1 = PG_GETARG_GSERIALIZED_P(0);
3135  geom2 = PG_GETARG_GSERIALIZED_P(1);
3136  tolerance = PG_GETARG_FLOAT8(2);
3137 
3138  lwgeom1 = lwgeom_from_gserialized(geom1);
3139  lwgeom2 = lwgeom_from_gserialized(geom2);
3140 
3141  lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3142  lwgeom_free(lwgeom1);
3143  lwgeom_free(lwgeom2);
3144  PG_FREE_IF_COPY(geom1, 0);
3145  PG_FREE_IF_COPY(geom2, 1);
3146 
3147  result = geometry_serialize(lwresult);
3148  lwgeom_free(lwresult);
3149 
3150  PG_RETURN_POINTER(result);
3151 }
3152 
3153 /*
3154  * ST_Split
3155  *
3156  * Split polygon by line, line by line, line by point.
3157  * Returns at most components as a collection.
3158  * First element of the collection is always the part which
3159  * remains after the cut, while the second element is the
3160  * part which has been cut out. We arbitrarily take the part
3161  * on the *right* of cut lines as the part which has been cut out.
3162  * For a line cut by a point the part which remains is the one
3163  * from start of the line to the cut point.
3164  *
3165  *
3166  * Author: Sandro Santilli <strk@kbt.io>
3167  *
3168  * Work done for Faunalia (http://www.faunalia.it) with fundings
3169  * from Regione Toscana - Sistema Informativo per il Governo
3170  * del Territorio e dell'Ambiente (RT-SIGTA).
3171  *
3172  * Thanks to the PostGIS community for sharing poly/line ideas [1]
3173  *
3174  * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3175  *
3176  */
3177 Datum ST_Split(PG_FUNCTION_ARGS);
3179 Datum ST_Split(PG_FUNCTION_ARGS)
3180 {
3181  GSERIALIZED *in, *blade_in, *out;
3182  LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
3183 
3184  in = PG_GETARG_GSERIALIZED_P(0);
3185  lwgeom_in = lwgeom_from_gserialized(in);
3186 
3187  blade_in = PG_GETARG_GSERIALIZED_P(1);
3188  lwblade_in = lwgeom_from_gserialized(blade_in);
3189 
3190  error_if_srid_mismatch(lwgeom_in->srid, lwblade_in->srid);
3191 
3192  lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3193  lwgeom_free(lwgeom_in);
3194  lwgeom_free(lwblade_in);
3195 
3196  if ( ! lwgeom_out )
3197  {
3198  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3199  PG_FREE_IF_COPY(blade_in, 1);
3200  PG_RETURN_NULL();
3201  }
3202 
3203  out = geometry_serialize(lwgeom_out);
3204  lwgeom_free(lwgeom_out);
3205  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3206  PG_FREE_IF_COPY(blade_in, 1);
3207 
3208  PG_RETURN_POINTER(out);
3209 }
3210 
3211 /**********************************************************************
3212  *
3213  * ST_SharedPaths
3214  *
3215  * Return the set of paths shared between two linear geometries,
3216  * and their direction (same or opposite).
3217  *
3218  * Developed by Sandro Santilli (strk@kbt.io) for Faunalia
3219  * (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3220  * Informativo per la Gestione del Territorio e dell' Ambiente
3221  * [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3222  * trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3223  * 0494241492)"
3224  *
3225  **********************************************************************/
3226 Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3228 Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3229 {
3230  GSERIALIZED *geom1, *geom2, *out;
3231  LWGEOM *g1, *g2, *lwgeom_out;
3232 
3233  geom1 = PG_GETARG_GSERIALIZED_P(0);
3234  geom2 = PG_GETARG_GSERIALIZED_P(1);
3235 
3236  g1 = lwgeom_from_gserialized(geom1);
3237  g2 = lwgeom_from_gserialized(geom2);
3238 
3239  lwgeom_out = lwgeom_sharedpaths(g1, g2);
3240  lwgeom_free(g1);
3241  lwgeom_free(g2);
3242 
3243  if ( ! lwgeom_out )
3244  {
3245  PG_FREE_IF_COPY(geom1, 0);
3246  PG_FREE_IF_COPY(geom2, 1);
3247  PG_RETURN_NULL();
3248  }
3249 
3250  out = geometry_serialize(lwgeom_out);
3251  lwgeom_free(lwgeom_out);
3252 
3253  PG_FREE_IF_COPY(geom1, 0);
3254  PG_FREE_IF_COPY(geom2, 1);
3255  PG_RETURN_POINTER(out);
3256 }
3257 
3258 
3259 /**********************************************************************
3260  *
3261  * ST_Node
3262  *
3263  * Fully node a set of lines using the least possible nodes while
3264  * preserving all of the input ones.
3265  *
3266  **********************************************************************/
3267 Datum ST_Node(PG_FUNCTION_ARGS);
3269 Datum ST_Node(PG_FUNCTION_ARGS)
3270 {
3271  GSERIALIZED *geom1, *out;
3272  LWGEOM *g1, *lwgeom_out;
3273 
3274  geom1 = PG_GETARG_GSERIALIZED_P(0);
3275 
3276  g1 = lwgeom_from_gserialized(geom1);
3277 
3278  lwgeom_out = lwgeom_node(g1);
3279  lwgeom_free(g1);
3280 
3281  if ( ! lwgeom_out )
3282  {
3283  PG_FREE_IF_COPY(geom1, 0);
3284  PG_RETURN_NULL();
3285  }
3286 
3287  out = geometry_serialize(lwgeom_out);
3288  lwgeom_free(lwgeom_out);
3289 
3290  PG_FREE_IF_COPY(geom1, 0);
3291  PG_RETURN_POINTER(out);
3292 }
3293 
3294 /******************************************
3295  *
3296  * ST_Voronoi
3297  *
3298  * Returns a Voronoi diagram constructed
3299  * from the points of the input geometry.
3300  *
3301  ******************************************/
3302 Datum ST_Voronoi(PG_FUNCTION_ARGS);
3304 Datum ST_Voronoi(PG_FUNCTION_ARGS)
3305 {
3306  GSERIALIZED* input;
3307  GSERIALIZED* clip;
3308  GSERIALIZED* result;
3309  LWGEOM* lwgeom_input;
3310  LWGEOM* lwgeom_result;
3311  double tolerance;
3312  GBOX clip_envelope;
3313  int custom_clip_envelope;
3314  int return_polygons;
3315 
3316  /* Return NULL on NULL geometry */
3317  if (PG_ARGISNULL(0))
3318  PG_RETURN_NULL();
3319 
3320  /* Read our tolerance value */
3321  if (PG_ARGISNULL(2))
3322  {
3323  lwpgerror("Tolerance must be a positive number.");
3324  PG_RETURN_NULL();
3325  }
3326 
3327  tolerance = PG_GETARG_FLOAT8(2);
3328 
3329  if (tolerance < 0)
3330  {
3331  lwpgerror("Tolerance must be a positive number.");
3332  PG_RETURN_NULL();
3333  }
3334 
3335  /* Are we returning lines or polygons? */
3336  if (PG_ARGISNULL(3))
3337  {
3338  lwpgerror("return_polygons must be true or false.");
3339  PG_RETURN_NULL();
3340  }
3341  return_polygons = PG_GETARG_BOOL(3);
3342 
3343  /* Read our clipping envelope, if applicable. */
3344  custom_clip_envelope = !PG_ARGISNULL(1);
3345  if (custom_clip_envelope) {
3346  clip = PG_GETARG_GSERIALIZED_P(1);
3347  if (!gserialized_get_gbox_p(clip, &clip_envelope))
3348  {
3349  lwpgerror("Could not determine envelope of clipping geometry.");
3350  PG_FREE_IF_COPY(clip, 1);
3351  PG_RETURN_NULL();
3352  }
3353  PG_FREE_IF_COPY(clip, 1);
3354  }
3355 
3356  /* Read our input geometry */
3357  input = PG_GETARG_GSERIALIZED_P(0);
3358 
3359  lwgeom_input = lwgeom_from_gserialized(input);
3360 
3361  if(!lwgeom_input)
3362  {
3363  lwpgerror("Could not read input geometry.");
3364  PG_FREE_IF_COPY(input, 0);
3365  PG_RETURN_NULL();
3366  }
3367 
3368  lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
3369  lwgeom_free(lwgeom_input);
3370 
3371  if (!lwgeom_result)
3372  {
3373  lwpgerror("Error computing Voronoi diagram.");
3374  PG_FREE_IF_COPY(input, 0);
3375  PG_RETURN_NULL();
3376  }
3377 
3378  result = geometry_serialize(lwgeom_result);
3379  lwgeom_free(lwgeom_result);
3380 
3381  PG_FREE_IF_COPY(input, 0);
3382  PG_RETURN_POINTER(result);
3383 }
3384 
3385 /******************************************
3386  *
3387  * ST_MinimumClearance
3388  *
3389  * Returns the minimum clearance of a geometry.
3390  *
3391  ******************************************/
3392 Datum ST_MinimumClearance(PG_FUNCTION_ARGS);
3394 Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
3395 {
3396 #if POSTGIS_GEOS_VERSION < 36
3397  lwpgerror("The GEOS version this PostGIS binary "
3398  "was compiled against (%d) doesn't support "
3399  "'ST_MinimumClearance' function (3.6.0+ required)",
3401  PG_RETURN_NULL();
3402 #else /* POSTGIS_GEOS_VERSION >= 36 */
3403  GSERIALIZED* input;
3404  GEOSGeometry* input_geos;
3405  int error;
3406  double result;
3407 
3408  initGEOS(lwpgnotice, lwgeom_geos_error);
3409 
3410  input = PG_GETARG_GSERIALIZED_P(0);
3411  input_geos = POSTGIS2GEOS(input);
3412  if (!input_geos)
3413  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3414 
3415  error = GEOSMinimumClearance(input_geos, &result);
3416  GEOSGeom_destroy(input_geos);
3417  if (error) HANDLE_GEOS_ERROR("Error computing minimum clearance");
3418 
3419  PG_FREE_IF_COPY(input, 0);
3420  PG_RETURN_FLOAT8(result);
3421 #endif
3422 }
3423 
3424 /******************************************
3425  *
3426  * ST_MinimumClearanceLine
3427  *
3428  * Returns the minimum clearance line of a geometry.
3429  *
3430  ******************************************/
3431 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS);
3433 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
3434 {
3435 #if POSTGIS_GEOS_VERSION < 36
3436  lwpgerror("The GEOS version this PostGIS binary "
3437  "was compiled against (%d) doesn't support "
3438  "'ST_MinimumClearanceLine' function (3.6.0+ required)",
3440  PG_RETURN_NULL();
3441 #else /* POSTGIS_GEOS_VERSION >= 36 */
3442  GSERIALIZED* input;
3443  GSERIALIZED* result;
3444  GEOSGeometry* input_geos;
3445  GEOSGeometry* result_geos;
3446  int srid;
3447 
3448  initGEOS(lwpgnotice, lwgeom_geos_error);
3449 
3450  input = PG_GETARG_GSERIALIZED_P(0);
3451  srid = gserialized_get_srid(input);
3452  input_geos = POSTGIS2GEOS(input);
3453  if (!input_geos)
3454  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3455 
3456  result_geos = GEOSMinimumClearanceLine(input_geos);
3457  GEOSGeom_destroy(input_geos);
3458  if (!result_geos)
3459  HANDLE_GEOS_ERROR("Error computing minimum clearance");
3460 
3461  GEOSSetSRID(result_geos, srid);
3462  result = GEOS2POSTGIS(result_geos, LW_FALSE);
3463  GEOSGeom_destroy(result_geos);
3464 
3465  PG_FREE_IF_COPY(input, 0);
3466  PG_RETURN_POINTER(result);
3467 #endif
3468 }
3469 
3470 /******************************************
3471  *
3472  * ST_OrientedEnvelope
3473  *
3474  ******************************************/
3475 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS);
3477 Datum ST_OrientedEnvelope(PG_FUNCTION_ARGS)
3478 {
3479 #if POSTGIS_GEOS_VERSION < 36
3480  lwpgerror("The GEOS version this PostGIS binary "
3481  "was compiled against (%d) doesn't support "
3482  "'ST_OrientedEnvelope' function (3.6.0+ required)",
3484  PG_RETURN_NULL();
3485 #else
3486  GSERIALIZED* input;
3487  GSERIALIZED* result;
3488  GEOSGeometry* input_geos;
3489  GEOSGeometry* result_geos;
3490  int srid;
3491 
3492  initGEOS(lwpgnotice, lwgeom_geos_error);
3493 
3494  input = PG_GETARG_GSERIALIZED_P(0);
3495  srid = gserialized_get_srid(input);
3496  input_geos = POSTGIS2GEOS(input);
3497  if (!input_geos)
3498  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3499 
3500  result_geos = GEOSMinimumRotatedRectangle(input_geos);
3501  GEOSGeom_destroy(input_geos);
3502  if (!result_geos)
3503  HANDLE_GEOS_ERROR("Error computing oriented envelope");
3504 
3505  GEOSSetSRID(result_geos, srid);
3506  result = GEOS2POSTGIS(result_geos, LW_FALSE);
3507  GEOSGeom_destroy(result_geos);
3508 
3509  PG_FREE_IF_COPY(input, 0);
3510  PG_RETURN_POINTER(result);
3511 #endif
3512 }
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:675
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.