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