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 = cstring2text(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  GEOSGeometry *g1, *g3;
815  GSERIALIZED *result;
816  int quadsegs = 8; /* the default */
817  int nargs;
818  enum
819  {
820  ENDCAP_ROUND = 1,
821  ENDCAP_FLAT = 2,
822  ENDCAP_SQUARE = 3
823  };
824  enum
825  {
826  JOIN_ROUND = 1,
827  JOIN_MITRE = 2,
828  JOIN_BEVEL = 3
829  };
830  static const double DEFAULT_MITRE_LIMIT = 5.0;
831  static const int DEFAULT_ENDCAP_STYLE = ENDCAP_ROUND;
832  static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
833 
834  double mitreLimit = DEFAULT_MITRE_LIMIT;
835  int endCapStyle = DEFAULT_ENDCAP_STYLE;
836  int joinStyle = DEFAULT_JOIN_STYLE;
837  char *param;
838  char *params = NULL;
839  LWGEOM *lwg;
840 
841  geom1 = PG_GETARG_GSERIALIZED_P(0);
842  size = PG_GETARG_FLOAT8(1);
843 
844  /* Empty.Buffer() == Empty[polygon] */
845  if ( gserialized_is_empty(geom1) )
846  {
848  gserialized_get_srid(geom1),
849  0, 0)); // buffer wouldn't give back z or m anyway
850  PG_RETURN_POINTER(geometry_serialize(lwg));
851  }
852 
853  nargs = PG_NARGS();
854 
855  initGEOS(lwpgnotice, lwgeom_geos_error);
856 
857  g1 = POSTGIS2GEOS(geom1);
858  if (!g1)
859  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
860 
861  if (nargs > 2)
862  {
863  /* We strdup `cause we're going to modify it */
864  params = pstrdup(PG_GETARG_CSTRING(2));
865 
866  POSTGIS_DEBUGF(3, "Params: %s", params);
867 
868  for (param=params; ; param=NULL)
869  {
870  char *key, *val;
871  param = strtok(param, " ");
872  if (!param) break;
873  POSTGIS_DEBUGF(3, "Param: %s", param);
874 
875  key = param;
876  val = strchr(key, '=');
877  if (!val || *(val + 1) == '\0')
878  {
879  lwpgerror("Missing value for buffer "
880  "parameter %s", key);
881  break;
882  }
883  *val = '\0';
884  ++val;
885 
886  POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
887 
888  if ( !strcmp(key, "endcap") )
889  {
890  /* Supported end cap styles:
891  * "round", "flat", "square"
892  */
893  if ( !strcmp(val, "round") )
894  {
895  endCapStyle = ENDCAP_ROUND;
896  }
897  else if ( !strcmp(val, "flat") ||
898  !strcmp(val, "butt") )
899  {
900  endCapStyle = ENDCAP_FLAT;
901  }
902  else if ( !strcmp(val, "square") )
903  {
904  endCapStyle = ENDCAP_SQUARE;
905  }
906  else
907  {
908  lwpgerror("Invalid buffer end cap "
909  "style: %s (accept: "
910  "'round', 'flat', 'butt' "
911  "or 'square'"
912  ")", val);
913  break;
914  }
915 
916  }
917  else if ( !strcmp(key, "join") )
918  {
919  if ( !strcmp(val, "round") )
920  {
921  joinStyle = JOIN_ROUND;
922  }
923  else if ( !strcmp(val, "mitre") ||
924  !strcmp(val, "miter") )
925  {
926  joinStyle = JOIN_MITRE;
927  }
928  else if ( !strcmp(val, "bevel") )
929  {
930  joinStyle = JOIN_BEVEL;
931  }
932  else
933  {
934  lwpgerror("Invalid buffer end cap "
935  "style: %s (accept: "
936  "'round', 'mitre', 'miter' "
937  " or 'bevel'"
938  ")", val);
939  break;
940  }
941  }
942  else if ( !strcmp(key, "mitre_limit") ||
943  !strcmp(key, "miter_limit") )
944  {
945  /* mitreLimit is a float */
946  mitreLimit = atof(val);
947  }
948  else if ( !strcmp(key, "quad_segs") )
949  {
950  /* quadrant segments is an int */
951  quadsegs = atoi(val);
952  }
953  else
954  {
955  lwpgerror("Invalid buffer parameter: %s (accept: "
956  "'endcap', 'join', 'mitre_limit', "
957  "'miter_limit and "
958  "'quad_segs')", key);
959  break;
960  }
961  }
962 
963  pfree(params); /* was pstrduped */
964 
965  POSTGIS_DEBUGF(3, "endCap:%d joinStyle:%d mitreLimit:%g",
966  endCapStyle, joinStyle, mitreLimit);
967 
968  }
969 
970  g3 = GEOSBufferWithStyle(g1, size, quadsegs, endCapStyle, joinStyle, mitreLimit);
971  GEOSGeom_destroy(g1);
972 
973  if (!g3) HANDLE_GEOS_ERROR("GEOSBuffer");
974 
975  POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3));
976 
977  GEOSSetSRID(g3, gserialized_get_srid(geom1));
978 
979  result = GEOS2POSTGIS(g3, gserialized_has_z(geom1));
980  GEOSGeom_destroy(g3);
981 
982  if (!result)
983  {
984  elog(ERROR,"GEOS buffer() threw an error (result postgis geometry formation)!");
985  PG_RETURN_NULL(); /* never get here */
986  }
987 
988  PG_FREE_IF_COPY(geom1, 0);
989  PG_RETURN_POINTER(result);
990 }
991 
992 
993 
994 /*
995 * Generate a field of random points within the area of a
996 * polygon or multipolygon. Throws an error for other geometry
997 * types.
998 */
999 Datum ST_GeneratePoints(PG_FUNCTION_ARGS);
1001 Datum ST_GeneratePoints(PG_FUNCTION_ARGS)
1002 {
1003  GSERIALIZED *gser_input;
1004  GSERIALIZED *gser_result;
1005  LWGEOM *lwgeom_input;
1006  LWGEOM *lwgeom_result;
1007  int32 npoints;
1008 
1009  gser_input = PG_GETARG_GSERIALIZED_P(0);
1010  npoints = DatumGetInt32(DirectFunctionCall1(numeric_int4, PG_GETARG_DATUM(1)));
1011 
1012  /* Smartasses get nothing back */
1013  if (npoints < 0)
1014  PG_RETURN_NULL();
1015 
1016  /* Types get checked in the code, we'll keep things small out there */
1017  lwgeom_input = lwgeom_from_gserialized(gser_input);
1018  lwgeom_result = (LWGEOM*)lwgeom_to_points(lwgeom_input, npoints);
1019  lwgeom_free(lwgeom_input);
1020  PG_FREE_IF_COPY(gser_input, 0);
1021 
1022  /* Return null as null */
1023  if (!lwgeom_result)
1024  PG_RETURN_NULL();
1025 
1026  /* Serialize and return */
1027  gser_result = gserialized_from_lwgeom(lwgeom_result, 0);
1028  lwgeom_free(lwgeom_result);
1029  PG_RETURN_POINTER(gser_result);
1030 }
1031 
1032 
1033 /*
1034 * Compute at offset curve to a line
1035 */
1036 Datum ST_OffsetCurve(PG_FUNCTION_ARGS);
1038 Datum ST_OffsetCurve(PG_FUNCTION_ARGS)
1039 {
1040  GSERIALIZED *gser_input;
1041  GSERIALIZED *gser_result;
1042  LWGEOM *lwgeom_input;
1043  LWGEOM *lwgeom_result;
1044  double size;
1045  int quadsegs = 8; /* the default */
1046  int nargs;
1047 
1048  enum
1049  {
1050  JOIN_ROUND = 1,
1051  JOIN_MITRE = 2,
1052  JOIN_BEVEL = 3
1053  };
1054 
1055  static const double DEFAULT_MITRE_LIMIT = 5.0;
1056  static const int DEFAULT_JOIN_STYLE = JOIN_ROUND;
1057  double mitreLimit = DEFAULT_MITRE_LIMIT;
1058  int joinStyle = DEFAULT_JOIN_STYLE;
1059  char *param = NULL;
1060  char *paramstr = NULL;
1061 
1062  /* Read SQL arguments */
1063  nargs = PG_NARGS();
1064  gser_input = PG_GETARG_GSERIALIZED_P(0);
1065  size = PG_GETARG_FLOAT8(1);
1066 
1067  /* Check for a useable type */
1068  if ( gserialized_get_type(gser_input) != LINETYPE )
1069  {
1070  lwpgerror("ST_OffsetCurve only works with LineStrings");
1071  PG_RETURN_NULL();
1072  }
1073 
1074  /*
1075  * For distance == 0, just return the input.
1076  * Note that due to a bug, GEOS 3.3.0 would return EMPTY.
1077  * See http://trac.osgeo.org/geos/ticket/454
1078  */
1079  if ( size == 0 )
1080  PG_RETURN_POINTER(gser_input);
1081 
1082  /* Read the lwgeom, check for errors */
1083  lwgeom_input = lwgeom_from_gserialized(gser_input);
1084  if ( ! lwgeom_input )
1085  lwpgerror("ST_OffsetCurve: lwgeom_from_gserialized returned NULL");
1086 
1087  /* For empty inputs, just echo them back */
1088  if ( lwgeom_is_empty(lwgeom_input) )
1089  PG_RETURN_POINTER(gser_input);
1090 
1091  /* Process the optional arguments */
1092  if ( nargs > 2 )
1093  {
1094  text *wkttext = PG_GETARG_TEXT_P(2);
1095  paramstr = text2cstring(wkttext);
1096 
1097  POSTGIS_DEBUGF(3, "paramstr: %s", paramstr);
1098 
1099  for ( param=paramstr; ; param=NULL )
1100  {
1101  char *key, *val;
1102  param = strtok(param, " ");
1103  if (!param) break;
1104  POSTGIS_DEBUGF(3, "Param: %s", param);
1105 
1106  key = param;
1107  val = strchr(key, '=');
1108  if (!val || *(val + 1) == '\0')
1109  {
1110  lwpgerror("ST_OffsetCurve: Missing value for buffer parameter %s", key);
1111  break;
1112  }
1113  *val = '\0';
1114  ++val;
1115 
1116  POSTGIS_DEBUGF(3, "Param: %s : %s", key, val);
1117 
1118  if ( !strcmp(key, "join") )
1119  {
1120  if ( !strcmp(val, "round") )
1121  {
1122  joinStyle = JOIN_ROUND;
1123  }
1124  else if ( !(strcmp(val, "mitre") && strcmp(val, "miter")) )
1125  {
1126  joinStyle = JOIN_MITRE;
1127  }
1128  else if ( ! strcmp(val, "bevel") )
1129  {
1130  joinStyle = JOIN_BEVEL;
1131  }
1132  else
1133  {
1134  lwpgerror("Invalid buffer end cap style: %s (accept: "
1135  "'round', 'mitre', 'miter' or 'bevel')", val);
1136  break;
1137  }
1138  }
1139  else if ( !strcmp(key, "mitre_limit") ||
1140  !strcmp(key, "miter_limit") )
1141  {
1142  /* mitreLimit is a float */
1143  mitreLimit = atof(val);
1144  }
1145  else if ( !strcmp(key, "quad_segs") )
1146  {
1147  /* quadrant segments is an int */
1148  quadsegs = atoi(val);
1149  }
1150  else
1151  {
1152  lwpgerror("Invalid buffer parameter: %s (accept: "
1153  "'join', 'mitre_limit', 'miter_limit and "
1154  "'quad_segs')", key);
1155  break;
1156  }
1157  }
1158  POSTGIS_DEBUGF(3, "joinStyle:%d mitreLimit:%g", joinStyle, mitreLimit);
1159  pfree(paramstr); /* alloc'ed in text2cstring */
1160  }
1161 
1162  lwgeom_result = lwgeom_offsetcurve(lwgeom_as_lwline(lwgeom_input), size, quadsegs, joinStyle, mitreLimit);
1163 
1164  if (!lwgeom_result)
1165  lwpgerror("ST_OffsetCurve: lwgeom_offsetcurve returned NULL");
1166 
1167  gser_result = gserialized_from_lwgeom(lwgeom_result, 0);
1168  lwgeom_free(lwgeom_input);
1169  lwgeom_free(lwgeom_result);
1170  PG_RETURN_POINTER(gser_result);
1171 }
1172 
1173 
1175 Datum geos_intersection(PG_FUNCTION_ARGS)
1176 {
1177  GSERIALIZED *geom1;
1178  GSERIALIZED *geom2;
1179  GSERIALIZED *result;
1180  LWGEOM *lwgeom1, *lwgeom2, *lwresult ;
1181 
1182  geom1 = PG_GETARG_GSERIALIZED_P(0);
1183  geom2 = PG_GETARG_GSERIALIZED_P(1);
1184 
1185  lwgeom1 = lwgeom_from_gserialized(geom1) ;
1186  lwgeom2 = lwgeom_from_gserialized(geom2) ;
1187 
1188  lwresult = lwgeom_intersection(lwgeom1, lwgeom2) ;
1189  result = geometry_serialize(lwresult) ;
1190 
1191  lwgeom_free(lwgeom1) ;
1192  lwgeom_free(lwgeom2) ;
1193  lwgeom_free(lwresult) ;
1194 
1195  PG_FREE_IF_COPY(geom1, 0);
1196  PG_FREE_IF_COPY(geom2, 1);
1197 
1198  PG_RETURN_POINTER(result);
1199 }
1200 
1207 Datum geos_difference(PG_FUNCTION_ARGS)
1208 {
1209  GSERIALIZED *geom1;
1210  GSERIALIZED *geom2;
1211  GSERIALIZED *result;
1212  LWGEOM *lwgeom1, *lwgeom2, *lwresult ;
1213 
1214  geom1 = PG_GETARG_GSERIALIZED_P(0);
1215  geom2 = PG_GETARG_GSERIALIZED_P(1);
1216 
1217  lwgeom1 = lwgeom_from_gserialized(geom1) ;
1218  lwgeom2 = lwgeom_from_gserialized(geom2) ;
1219 
1220  lwresult = lwgeom_difference(lwgeom1, lwgeom2) ;
1221  result = geometry_serialize(lwresult) ;
1222 
1223  lwgeom_free(lwgeom1) ;
1224  lwgeom_free(lwgeom2) ;
1225  lwgeom_free(lwresult) ;
1226 
1227  PG_FREE_IF_COPY(geom1, 0);
1228  PG_FREE_IF_COPY(geom2, 1);
1229 
1230  PG_RETURN_POINTER(result);
1231 }
1232 
1233 
1238 Datum pointonsurface(PG_FUNCTION_ARGS)
1239 {
1240  GSERIALIZED *geom;
1241  GEOSGeometry *g1, *g3;
1242  GSERIALIZED *result;
1243 
1244  geom = PG_GETARG_GSERIALIZED_P(0);
1245 
1246  /* Empty.PointOnSurface == Point Empty */
1247  if ( gserialized_is_empty(geom) )
1248  {
1250  gserialized_get_srid(geom),
1251  gserialized_has_z(geom),
1252  gserialized_has_m(geom));
1253  result = geometry_serialize(lwpoint_as_lwgeom(lwp));
1254  lwpoint_free(lwp);
1255  PG_RETURN_POINTER(result);
1256  }
1257 
1258  initGEOS(lwpgnotice, lwgeom_geos_error);
1259 
1260  g1 = POSTGIS2GEOS(geom);
1261 
1262  if (!g1)
1263  {
1264  /* Why is this a WARNING rather than an error ? */
1265  /* TODO: use HANDLE_GEOS_ERROR instead */
1266  elog(WARNING, "GEOSPointOnSurface(): %s", lwgeom_geos_errmsg);
1267  PG_RETURN_NULL();
1268  }
1269 
1270  g3 = GEOSPointOnSurface(g1);
1271 
1272  if (!g3)
1273  {
1274  GEOSGeom_destroy(g1);
1275  HANDLE_GEOS_ERROR("GEOSPointOnSurface");
1276  }
1277 
1278  POSTGIS_DEBUGF(3, "result: %s", GEOSGeomToWKT(g3) ) ;
1279 
1280  GEOSSetSRID(g3, gserialized_get_srid(geom));
1281 
1282  result = GEOS2POSTGIS(g3, gserialized_has_z(geom));
1283 
1284  if (!result)
1285  {
1286  GEOSGeom_destroy(g1);
1287  GEOSGeom_destroy(g3);
1288  elog(ERROR,"GEOS pointonsurface() threw an error (result postgis geometry formation)!");
1289  PG_RETURN_NULL(); /* never get here */
1290  }
1291 
1292  GEOSGeom_destroy(g1);
1293  GEOSGeom_destroy(g3);
1294 
1295  PG_FREE_IF_COPY(geom, 0);
1296 
1297  PG_RETURN_POINTER(result);
1298 }
1299 
1301 Datum centroid(PG_FUNCTION_ARGS)
1302 {
1303  GSERIALIZED *geom, *result;
1304  GEOSGeometry *geosgeom, *geosresult;
1305  LWGEOM *igeom = NULL, *linear_geom = NULL;
1306  int32 perQuad= 16;
1307  int type = 0;
1308  geom = PG_GETARG_GSERIALIZED_P(0);
1309 
1310  /* Empty.Centroid() == Point Empty */
1311  if ( gserialized_is_empty(geom) )
1312  {
1314  gserialized_get_srid(geom),
1315  gserialized_has_z(geom),
1316  gserialized_has_m(geom));
1317  result = geometry_serialize(lwpoint_as_lwgeom(lwp));
1318  lwpoint_free(lwp);
1319  PG_RETURN_POINTER(result);
1320  }
1321 
1322  type = gserialized_get_type(geom) ;
1323  /* Converting curve geometry to linestring if necessary*/
1324  if(type == CIRCSTRINGTYPE || type == COMPOUNDTYPE )
1325  {/* curve geometry?*/
1326  igeom = lwgeom_from_gserialized(geom);
1327  PG_FREE_IF_COPY(geom, 0); /*free memory, we already have a lwgeom geometry copy*/
1328  linear_geom = lwgeom_stroke(igeom, perQuad);
1329  lwgeom_free(igeom);
1330  if (!linear_geom) PG_RETURN_NULL();
1331 
1332  geom = geometry_serialize(linear_geom);
1333  lwgeom_free(linear_geom);
1334  }
1335 
1336  initGEOS(lwpgnotice, lwgeom_geos_error);
1337 
1338  geosgeom = POSTGIS2GEOS(geom);
1339 
1340  if (!geosgeom)
1341  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1342 
1343  geosresult = GEOSGetCentroid(geosgeom);
1344 
1345  if (!geosresult)
1346  {
1347  GEOSGeom_destroy(geosgeom);
1348  HANDLE_GEOS_ERROR("GEOSGetCentroid");
1349  }
1350 
1351  GEOSSetSRID(geosresult, gserialized_get_srid(geom));
1352 
1353  result = GEOS2POSTGIS(geosresult, gserialized_has_z(geom));
1354 
1355  if (!result)
1356  {
1357  GEOSGeom_destroy(geosgeom);
1358  GEOSGeom_destroy(geosresult);
1359  elog(ERROR,"Error in GEOS-PGIS conversion");
1360  PG_RETURN_NULL();
1361  }
1362  GEOSGeom_destroy(geosgeom);
1363  GEOSGeom_destroy(geosresult);
1364 
1365  PG_FREE_IF_COPY(geom, 0);
1366 
1367  PG_RETURN_POINTER(result);
1368 }
1369 
1370 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS);
1372 Datum ST_ClipByBox2d(PG_FUNCTION_ARGS)
1373 {
1374 #if POSTGIS_GEOS_VERSION < 35
1375 
1376  lwpgerror("The GEOS version this PostGIS binary "
1377  "was compiled against (%d) doesn't support "
1378  "'GEOSClipByRect' function (3.5.0+ required)",
1380  PG_RETURN_NULL();
1381 
1382 #else /* POSTGIS_GEOS_VERSION >= 35 */
1383 
1384  GSERIALIZED *geom1;
1385  GSERIALIZED *result;
1386  LWGEOM *lwgeom1, *lwresult ;
1387  const GBOX *bbox1;
1388  GBOX *bbox2;
1389 
1390  geom1 = PG_GETARG_GSERIALIZED_P(0);
1391  lwgeom1 = lwgeom_from_gserialized(geom1) ;
1392 
1393  bbox1 = lwgeom_get_bbox(lwgeom1);
1394  if ( ! bbox1 )
1395  {
1396  /* empty clips to empty, no matter rect */
1397  lwgeom_free(lwgeom1);
1398  PG_RETURN_POINTER(geom1);
1399  }
1400 
1401  /* WARNING: this is really a BOX2DF, use only xmin and ymin fields */
1402  bbox2 = (GBOX *)PG_GETARG_POINTER(1);
1403  bbox2->flags = 0;
1404 
1405  /* If bbox1 outside of bbox2, return empty */
1406  if ( ! gbox_overlaps_2d(bbox1, bbox2) )
1407  {
1408  lwresult = lwgeom_construct_empty(lwgeom1->type, lwgeom1->srid, 0, 0);
1409  lwgeom_free(lwgeom1);
1410  PG_FREE_IF_COPY(geom1, 0);
1411  result = geometry_serialize(lwresult) ;
1412  lwgeom_free(lwresult) ;
1413  PG_RETURN_POINTER(result);
1414  }
1415 
1416  /* if bbox1 is covered by bbox2, return lwgeom1 */
1417  if ( gbox_contains_2d(bbox2, bbox1) )
1418  {
1419  lwgeom_free(lwgeom1);
1420  PG_RETURN_POINTER(geom1);
1421  }
1422 
1423  lwresult = lwgeom_clip_by_rect(lwgeom1, bbox2->xmin, bbox2->ymin,
1424  bbox2->xmax, bbox2->ymax);
1425 
1426  lwgeom_free(lwgeom1);
1427  PG_FREE_IF_COPY(geom1, 0);
1428 
1429  if (!lwresult) PG_RETURN_NULL();
1430 
1431  result = geometry_serialize(lwresult) ;
1432  lwgeom_free(lwresult) ;
1433  PG_RETURN_POINTER(result);
1434 
1435 #endif /* POSTGIS_GEOS_VERSION >= 35 */
1436 }
1437 
1438 
1439 
1440 /*---------------------------------------------*/
1441 
1442 
1450 {
1451  int t1 = gserialized_get_type(g1);
1452  int t2 = gserialized_get_type(g2);
1453 
1454  char *hintmsg;
1455  char *hintwkt;
1456  size_t hintsz;
1457  LWGEOM *lwgeom;
1458 
1459  if ( t1 == COLLECTIONTYPE)
1460  {
1461  lwgeom = lwgeom_from_gserialized(g1);
1462  hintwkt = lwgeom_to_wkt(lwgeom, WKT_SFSQL, DBL_DIG, &hintsz);
1463  lwgeom_free(lwgeom);
1464  hintmsg = lwmessage_truncate(hintwkt, 0, hintsz-1, 80, 1);
1465  ereport(ERROR,
1466  (errmsg("Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported."),
1467  errhint("Change argument 1: '%s'", hintmsg))
1468  );
1469  pfree(hintwkt);
1470  pfree(hintmsg);
1471  }
1472  else if (t2 == COLLECTIONTYPE)
1473  {
1474  lwgeom = lwgeom_from_gserialized(g2);
1475  hintwkt = lwgeom_to_wkt(lwgeom, WKT_SFSQL, DBL_DIG, &hintsz);
1476  hintmsg = lwmessage_truncate(hintwkt, 0, hintsz-1, 80, 1);
1477  lwgeom_free(lwgeom);
1478  ereport(ERROR,
1479  (errmsg("Relate Operation called with a LWGEOMCOLLECTION type. This is unsupported."),
1480  errhint("Change argument 2: '%s'", hintmsg))
1481  );
1482  pfree(hintwkt);
1483  pfree(hintmsg);
1484  }
1485 }
1486 
1488 Datum isvalid(PG_FUNCTION_ARGS)
1489 {
1490  GSERIALIZED *geom1;
1491  LWGEOM *lwgeom;
1492  bool result;
1493  GEOSGeom g1;
1494 
1495  geom1 = PG_GETARG_GSERIALIZED_P(0);
1496 
1497  /* Empty.IsValid() == TRUE */
1498  if ( gserialized_is_empty(geom1) )
1499  PG_RETURN_BOOL(true);
1500 
1501  initGEOS(lwpgnotice, lwgeom_geos_error);
1502 
1503  lwgeom = lwgeom_from_gserialized(geom1);
1504  if ( ! lwgeom )
1505  {
1506  lwpgerror("unable to deserialize input");
1507  }
1508  g1 = LWGEOM2GEOS(lwgeom, 0);
1509  lwgeom_free(lwgeom);
1510 
1511  if ( ! g1 )
1512  {
1513  /* should we drop the following
1514  * notice now that we have ST_isValidReason ?
1515  */
1516  lwpgnotice("%s", lwgeom_geos_errmsg);
1517  PG_RETURN_BOOL(false);
1518  }
1519 
1520  result = GEOSisValid(g1);
1521  GEOSGeom_destroy(g1);
1522 
1523  if (result == 2)
1524  {
1525  elog(ERROR,"GEOS isvalid() threw an error!");
1526  PG_RETURN_NULL(); /*never get here */
1527  }
1528 
1529  PG_FREE_IF_COPY(geom1, 0);
1530  PG_RETURN_BOOL(result);
1531 }
1532 
1533 /*
1534 ** IsValidReason is only available in the GEOS
1535 ** C API > version 3.0
1536 */
1538 Datum isvalidreason(PG_FUNCTION_ARGS)
1539 {
1540  GSERIALIZED *geom = NULL;
1541  char *reason_str = NULL;
1542  text *result = NULL;
1543  const GEOSGeometry *g1 = NULL;
1544 
1545  geom = PG_GETARG_GSERIALIZED_P(0);
1546 
1547  initGEOS(lwpgnotice, lwgeom_geos_error);
1548 
1549  g1 = POSTGIS2GEOS(geom);
1550  if ( g1 )
1551  {
1552  reason_str = GEOSisValidReason(g1);
1553  GEOSGeom_destroy((GEOSGeometry *)g1);
1554  if (!reason_str) HANDLE_GEOS_ERROR("GEOSisValidReason");
1555  result = cstring2text(reason_str);
1556  GEOSFree(reason_str);
1557  }
1558  else
1559  {
1560  result = cstring2text(lwgeom_geos_errmsg);
1561  }
1562 
1563  PG_FREE_IF_COPY(geom, 0);
1564  PG_RETURN_POINTER(result);
1565 }
1566 
1567 /*
1568 ** IsValidDetail is only available in the GEOS
1569 ** C API >= version 3.3
1570 */
1572 Datum isvaliddetail(PG_FUNCTION_ARGS)
1573 {
1574  GSERIALIZED *geom = NULL;
1575  const GEOSGeometry *g1 = NULL;
1576  char *values[3]; /* valid bool, reason text, location geometry */
1577  char *geos_reason = NULL;
1578  char *reason = NULL;
1579  GEOSGeometry *geos_location = NULL;
1580  LWGEOM *location = NULL;
1581  char valid = 0;
1582  HeapTupleHeader result;
1583  TupleDesc tupdesc;
1584  HeapTuple tuple;
1585  AttInMetadata *attinmeta;
1586  int flags = 0;
1587 
1588  /*
1589  * Build a tuple description for a
1590  * valid_detail tuple
1591  */
1592  tupdesc = RelationNameGetTupleDesc("valid_detail");
1593  if ( ! tupdesc )
1594  {
1595  lwpgerror("TYPE valid_detail not found");
1596  PG_RETURN_NULL();
1597  }
1598 
1599  /*
1600  * generate attribute metadata needed later to produce
1601  * tuples from raw C strings
1602  */
1603  attinmeta = TupleDescGetAttInMetadata(tupdesc);
1604 
1605  geom = PG_GETARG_GSERIALIZED_P(0);
1606 
1607  if ( PG_NARGS() > 1 && ! PG_ARGISNULL(1) )
1608  {
1609  flags = PG_GETARG_INT32(1);
1610  }
1611 
1612  initGEOS(lwpgnotice, lwgeom_geos_error);
1613 
1614  g1 = POSTGIS2GEOS(geom);
1615 
1616  if ( g1 )
1617  {
1618  valid = GEOSisValidDetail(g1, flags,
1619  &geos_reason, &geos_location);
1620  GEOSGeom_destroy((GEOSGeometry *)g1);
1621  if ( geos_reason )
1622  {
1623  reason = pstrdup(geos_reason);
1624  GEOSFree(geos_reason);
1625  }
1626  if ( geos_location )
1627  {
1628  location = GEOS2LWGEOM(geos_location, GEOSHasZ(geos_location));
1629  GEOSGeom_destroy(geos_location);
1630  }
1631 
1632  if (valid == 2)
1633  {
1634  /* NOTE: should only happen on OOM or similar */
1635  lwpgerror("GEOS isvaliddetail() threw an exception!");
1636  PG_RETURN_NULL(); /* never gets here */
1637  }
1638  }
1639  else
1640  {
1641  /* TODO: check lwgeom_geos_errmsg for validity error */
1642  reason = pstrdup(lwgeom_geos_errmsg);
1643  }
1644 
1645  /* the boolean validity */
1646  values[0] = valid ? "t" : "f";
1647 
1648  /* the reason */
1649  values[1] = reason;
1650 
1651  /* the location */
1652  values[2] = location ? lwgeom_to_hexwkb(location, WKB_EXTENDED, 0) : 0;
1653 
1654  tuple = BuildTupleFromCStrings(attinmeta, values);
1655  result = (HeapTupleHeader) palloc(tuple->t_len);
1656  memcpy(result, tuple->t_data, tuple->t_len);
1657  heap_freetuple(tuple);
1658 
1659  PG_RETURN_HEAPTUPLEHEADER(result);
1660 }
1661 
1670 Datum overlaps(PG_FUNCTION_ARGS)
1671 {
1672  GSERIALIZED *geom1;
1673  GSERIALIZED *geom2;
1674  GEOSGeometry *g1, *g2;
1675  bool result;
1676  GBOX box1, box2;
1677 
1678  geom1 = PG_GETARG_GSERIALIZED_P(0);
1679  geom2 = PG_GETARG_GSERIALIZED_P(1);
1680 
1681  errorIfGeometryCollection(geom1,geom2);
1683 
1684  /* A.Overlaps(Empty) == FALSE */
1685  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1686  PG_RETURN_BOOL(false);
1687 
1688  /*
1689  * short-circuit 1: if geom2 bounding box does not overlap
1690  * geom1 bounding box we can prematurely return FALSE.
1691  */
1692  if ( gserialized_get_gbox_p(geom1, &box1) &&
1693  gserialized_get_gbox_p(geom2, &box2) )
1694  {
1695  if ( ! gbox_overlaps_2d(&box1, &box2) )
1696  {
1697  PG_RETURN_BOOL(false);
1698  }
1699  }
1700 
1701  initGEOS(lwpgnotice, lwgeom_geos_error);
1702 
1703  g1 = POSTGIS2GEOS(geom1);
1704  if (!g1)
1705  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
1706 
1707  g2 = POSTGIS2GEOS(geom2);
1708 
1709  if (!g2)
1710  {
1711  GEOSGeom_destroy(g1);
1712  HANDLE_GEOS_ERROR("Second argument geometry could not be converted to GEOS");
1713  }
1714 
1715  result = GEOSOverlaps(g1,g2);
1716 
1717  GEOSGeom_destroy(g1);
1718  GEOSGeom_destroy(g2);
1719  if (result == 2) HANDLE_GEOS_ERROR("GEOSOverlaps");
1720 
1721  PG_FREE_IF_COPY(geom1, 0);
1722  PG_FREE_IF_COPY(geom2, 1);
1723 
1724  PG_RETURN_BOOL(result);
1725 }
1726 
1727 
1729 Datum contains(PG_FUNCTION_ARGS)
1730 {
1731  GSERIALIZED *geom1;
1732  GSERIALIZED *geom2;
1733  GEOSGeometry *g1, *g2;
1734  GBOX box1, box2;
1735  int result;
1736  PrepGeomCache *prep_cache;
1737 
1738  geom1 = PG_GETARG_GSERIALIZED_P(0);
1739  geom2 = PG_GETARG_GSERIALIZED_P(1);
1740 
1741  errorIfGeometryCollection(geom1,geom2);
1743 
1744  /* A.Contains(Empty) == FALSE */
1745  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1746  PG_RETURN_BOOL(false);
1747 
1748  POSTGIS_DEBUG(3, "contains called.");
1749 
1750  /*
1751  ** short-circuit 1: if geom2 bounding box is not completely inside
1752  ** geom1 bounding box we can prematurely return FALSE.
1753  */
1754  if ( gserialized_get_gbox_p(geom1, &box1) &&
1755  gserialized_get_gbox_p(geom2, &box2) )
1756  {
1757  if ( ! gbox_contains_2d(&box1, &box2) )
1758  {
1759  PG_RETURN_BOOL(false);
1760  }
1761  }
1762 
1763  /*
1764  ** short-circuit 2: if geom2 is a point and geom1 is a polygon
1765  ** call the point-in-polygon function.
1766  */
1767  if (is_poly(geom1) && is_point(geom2))
1768  {
1769  GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
1770  GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
1771  RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
1772  int retval;
1773 
1774  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1775  if (gserialized_get_type(gpoint) == POINTTYPE)
1776  {
1777  LWGEOM* point = lwgeom_from_gserialized(gpoint);
1778  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1779  lwgeom_free(point);
1780 
1781  retval = (pip_result == 1); /* completely inside */
1782  }
1783  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
1784  {
1786  uint32_t i;
1787  int found_completely_inside = LW_FALSE;
1788 
1789  retval = LW_TRUE;
1790  for (i = 0; i < mpoint->ngeoms; i++)
1791  {
1792  /* We need to find at least one point that's completely inside the
1793  * polygons (pip_result == 1). As long as we have one point that's
1794  * completely inside, we can have as many as we want on the boundary
1795  * itself. (pip_result == 0)
1796  */
1797  int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
1798  if (pip_result == 1)
1799  found_completely_inside = LW_TRUE;
1800 
1801  if (pip_result == -1) /* completely outside */
1802  {
1803  retval = LW_FALSE;
1804  break;
1805  }
1806  }
1807 
1808  retval = retval && found_completely_inside;
1809  lwmpoint_free(mpoint);
1810  }
1811  else
1812  {
1813  /* Never get here */
1814  elog(ERROR,"Type isn't point or multipoint!");
1815  PG_RETURN_NULL();
1816  }
1817 
1818  PG_FREE_IF_COPY(geom1, 0);
1819  PG_FREE_IF_COPY(geom2, 1);
1820  PG_RETURN_BOOL(retval);
1821  }
1822  else
1823  {
1824  POSTGIS_DEBUGF(3, "Contains: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
1825  }
1826 
1827  initGEOS(lwpgnotice, lwgeom_geos_error);
1828 
1829  prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
1830 
1831  if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
1832  {
1833  g1 = POSTGIS2GEOS(geom2);
1834  if (!g1)
1835  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
1836 
1837  POSTGIS_DEBUG(4, "containsPrepared: cache is live, running preparedcontains");
1838  result = GEOSPreparedContains( prep_cache->prepared_geom, g1);
1839  GEOSGeom_destroy(g1);
1840  }
1841  else
1842  {
1843  g1 = POSTGIS2GEOS(geom1);
1844  if (!g1)
1846  "First argument geometry could not be converted to "
1847  "GEOS");
1848  g2 = POSTGIS2GEOS(geom2);
1849  if (!g2)
1850  {
1852  "Second argument geometry could not be converted "
1853  "to GEOS");
1854  GEOSGeom_destroy(g1);
1855  }
1856  POSTGIS_DEBUG(4, "containsPrepared: cache is not ready, running standard contains");
1857  result = GEOSContains( g1, g2);
1858  GEOSGeom_destroy(g1);
1859  GEOSGeom_destroy(g2);
1860  }
1861 
1862  if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1863 
1864  PG_FREE_IF_COPY(geom1, 0);
1865  PG_FREE_IF_COPY(geom2, 1);
1866 
1867  PG_RETURN_BOOL(result);
1868 
1869 }
1870 
1872 Datum containsproperly(PG_FUNCTION_ARGS)
1873 {
1874  GSERIALIZED * geom1;
1875  GSERIALIZED * geom2;
1876  bool result;
1877  GBOX box1, box2;
1878  PrepGeomCache * prep_cache;
1879 
1880  geom1 = PG_GETARG_GSERIALIZED_P(0);
1881  geom2 = PG_GETARG_GSERIALIZED_P(1);
1882 
1883  errorIfGeometryCollection(geom1,geom2);
1885 
1886  /* A.ContainsProperly(Empty) == FALSE */
1887  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1888  PG_RETURN_BOOL(false);
1889 
1890  /*
1891  * short-circuit: if geom2 bounding box is not completely inside
1892  * geom1 bounding box we can prematurely return FALSE.
1893  */
1894  if ( gserialized_get_gbox_p(geom1, &box1) &&
1895  gserialized_get_gbox_p(geom2, &box2) )
1896  {
1897  if ( ! gbox_contains_2d(&box1, &box2) )
1898  PG_RETURN_BOOL(false);
1899  }
1900 
1901  initGEOS(lwpgnotice, lwgeom_geos_error);
1902 
1903  prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
1904 
1905  if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
1906  {
1907  GEOSGeometry *g = POSTGIS2GEOS(geom2);
1908  if (!g)
1910  "First argument geometry could not be converted to "
1911  "GEOS");
1912  result = GEOSPreparedContainsProperly( prep_cache->prepared_geom, g);
1913  GEOSGeom_destroy(g);
1914  }
1915  else
1916  {
1917  GEOSGeometry *g2;
1918  GEOSGeometry *g1;
1919 
1920  g1 = POSTGIS2GEOS(geom1);
1921  if (!g1)
1923  "First argument geometry could not be converted to "
1924  "GEOS");
1925  g2 = POSTGIS2GEOS(geom2);
1926  if (!g2)
1927  {
1928  GEOSGeom_destroy(g1);
1930  "Second argument geometry could not be converted "
1931  "to GEOS");
1932  }
1933  result = GEOSRelatePattern( g1, g2, "T**FF*FF*" );
1934  GEOSGeom_destroy(g1);
1935  GEOSGeom_destroy(g2);
1936  }
1937 
1938  if (result == 2) HANDLE_GEOS_ERROR("GEOSContains");
1939 
1940  PG_FREE_IF_COPY(geom1, 0);
1941  PG_FREE_IF_COPY(geom2, 1);
1942 
1943  PG_RETURN_BOOL(result);
1944 }
1945 
1946 /*
1947  * Described at
1948  * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
1949  */
1951 Datum covers(PG_FUNCTION_ARGS)
1952 {
1953  GSERIALIZED *geom1;
1954  GSERIALIZED *geom2;
1955  int result;
1956  GBOX box1, box2;
1957  PrepGeomCache *prep_cache;
1958 
1959  geom1 = PG_GETARG_GSERIALIZED_P(0);
1960  geom2 = PG_GETARG_GSERIALIZED_P(1);
1961 
1962  /* A.Covers(Empty) == FALSE */
1963  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
1964  PG_RETURN_BOOL(false);
1965 
1966  errorIfGeometryCollection(geom1,geom2);
1968 
1969  /*
1970  * short-circuit 1: if geom2 bounding box is not completely inside
1971  * geom1 bounding box we can prematurely return FALSE.
1972  */
1973  if ( gserialized_get_gbox_p(geom1, &box1) &&
1974  gserialized_get_gbox_p(geom2, &box2) )
1975  {
1976  if ( ! gbox_contains_2d(&box1, &box2) )
1977  {
1978  PG_RETURN_BOOL(false);
1979  }
1980  }
1981  /*
1982  * short-circuit 2: if geom2 is a point and geom1 is a polygon
1983  * call the point-in-polygon function.
1984  */
1985  if (is_poly(geom1) && is_point(geom2))
1986  {
1987  GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
1988  GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
1989  RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
1990  int retval;
1991 
1992  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
1993  if (gserialized_get_type(gpoint) == POINTTYPE)
1994  {
1995  LWGEOM* point = lwgeom_from_gserialized(gpoint);
1996  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
1997  lwgeom_free(point);
1998 
1999  retval = (pip_result != -1); /* not outside */
2000  }
2001  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2002  {
2004  uint32_t i;
2005 
2006  retval = LW_TRUE;
2007  for (i = 0; i < mpoint->ngeoms; i++)
2008  {
2009  int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2010  if (pip_result == -1)
2011  {
2012  retval = LW_FALSE;
2013  break;
2014  }
2015  }
2016 
2017  lwmpoint_free(mpoint);
2018  }
2019  else
2020  {
2021  /* Never get here */
2022  elog(ERROR,"Type isn't point or multipoint!");
2023  PG_RETURN_NULL();
2024  }
2025 
2026  PG_FREE_IF_COPY(geom1, 0);
2027  PG_FREE_IF_COPY(geom2, 1);
2028  PG_RETURN_BOOL(retval);
2029  }
2030  else
2031  {
2032  POSTGIS_DEBUGF(3, "Covers: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2033  }
2034 
2035  initGEOS(lwpgnotice, lwgeom_geos_error);
2036 
2037  prep_cache = GetPrepGeomCache( fcinfo, geom1, 0 );
2038 
2039  if ( prep_cache && prep_cache->prepared_geom && prep_cache->argnum == 1 )
2040  {
2041  GEOSGeometry *g1 = POSTGIS2GEOS(geom2);
2042  if (!g1)
2044  "First argument geometry could not be converted to "
2045  "GEOS");
2046  result = GEOSPreparedCovers( prep_cache->prepared_geom, g1);
2047  GEOSGeom_destroy(g1);
2048  }
2049  else
2050  {
2051  GEOSGeometry *g1;
2052  GEOSGeometry *g2;
2053 
2054  g1 = POSTGIS2GEOS(geom1);
2055  if (!g1)
2057  "First argument geometry could not be converted to "
2058  "GEOS");
2059  g2 = POSTGIS2GEOS(geom2);
2060  if (!g2)
2061  {
2062  GEOSGeom_destroy(g1);
2064  "Second argument geometry could not be converted "
2065  "to GEOS");
2066  }
2067  result = GEOSRelatePattern( g1, g2, "******FF*" );
2068  GEOSGeom_destroy(g1);
2069  GEOSGeom_destroy(g2);
2070  }
2071 
2072  if (result == 2) HANDLE_GEOS_ERROR("GEOSCovers");
2073 
2074  PG_FREE_IF_COPY(geom1, 0);
2075  PG_FREE_IF_COPY(geom2, 1);
2076 
2077  PG_RETURN_BOOL(result);
2078 
2079 }
2080 
2081 
2089 /*
2090  * Described at:
2091  * http://lin-ear-th-inking.blogspot.com/2007/06/subtleties-of-ogc-covers-spatial.html
2092  */
2094 Datum coveredby(PG_FUNCTION_ARGS)
2095 {
2096  GSERIALIZED *geom1;
2097  GSERIALIZED *geom2;
2098  GEOSGeometry *g1, *g2;
2099  int result;
2100  GBOX box1, box2;
2101  char *patt = "**F**F***";
2102 
2103  geom1 = PG_GETARG_GSERIALIZED_P(0);
2104  geom2 = PG_GETARG_GSERIALIZED_P(1);
2105 
2106  errorIfGeometryCollection(geom1,geom2);
2108 
2109  /* A.CoveredBy(Empty) == FALSE */
2110  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2111  PG_RETURN_BOOL(false);
2112 
2113  /*
2114  * short-circuit 1: if geom1 bounding box is not completely inside
2115  * geom2 bounding box we can prematurely return FALSE.
2116  */
2117  if ( gserialized_get_gbox_p(geom1, &box1) &&
2118  gserialized_get_gbox_p(geom2, &box2) )
2119  {
2120  if ( ! gbox_contains_2d(&box2, &box1) )
2121  {
2122  PG_RETURN_BOOL(false);
2123  }
2124 
2125  POSTGIS_DEBUG(3, "bounding box short-circuit missed.");
2126  }
2127  /*
2128  * short-circuit 2: if geom1 is a point and geom2 is a polygon
2129  * call the point-in-polygon function.
2130  */
2131  if (is_point(geom1) && is_poly(geom2))
2132  {
2133  GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
2134  GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
2135  RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
2136  int retval;
2137 
2138  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2139  if (gserialized_get_type(gpoint) == POINTTYPE)
2140  {
2141  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2142  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2143  lwgeom_free(point);
2144 
2145  retval = (pip_result != -1); /* not outside */
2146  }
2147  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2148  {
2150  uint32_t i;
2151 
2152  retval = LW_TRUE;
2153  for (i = 0; i < mpoint->ngeoms; i++)
2154  {
2155  int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2156  if (pip_result == -1)
2157  {
2158  retval = LW_FALSE;
2159  break;
2160  }
2161  }
2162 
2163  lwmpoint_free(mpoint);
2164  }
2165  else
2166  {
2167  /* Never get here */
2168  elog(ERROR,"Type isn't point or multipoint!");
2169  PG_RETURN_NULL();
2170  }
2171 
2172  PG_FREE_IF_COPY(geom1, 0);
2173  PG_FREE_IF_COPY(geom2, 1);
2174  PG_RETURN_BOOL(retval);
2175  }
2176  else
2177  {
2178  POSTGIS_DEBUGF(3, "CoveredBy: type1: %d, type2: %d", gserialized_get_type(geom1), gserialized_get_type(geom2));
2179  }
2180 
2181  initGEOS(lwpgnotice, lwgeom_geos_error);
2182 
2183  g1 = POSTGIS2GEOS(geom1);
2184 
2185  if (!g1)
2186  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2187 
2188  g2 = POSTGIS2GEOS(geom2);
2189 
2190  if (!g2)
2191  {
2192  GEOSGeom_destroy(g1);
2194  "Second argument geometry could not be converted to GEOS");
2195  }
2196 
2197  result = GEOSRelatePattern(g1,g2,patt);
2198 
2199  GEOSGeom_destroy(g1);
2200  GEOSGeom_destroy(g2);
2201 
2202  if (result == 2) HANDLE_GEOS_ERROR("GEOSCoveredBy");
2203 
2204  PG_FREE_IF_COPY(geom1, 0);
2205  PG_FREE_IF_COPY(geom2, 1);
2206 
2207  PG_RETURN_BOOL(result);
2208 }
2209 
2211 Datum crosses(PG_FUNCTION_ARGS)
2212 {
2213  GSERIALIZED *geom1;
2214  GSERIALIZED *geom2;
2215  GEOSGeometry *g1, *g2;
2216  int result;
2217  GBOX box1, box2;
2218 
2219  geom1 = PG_GETARG_GSERIALIZED_P(0);
2220  geom2 = PG_GETARG_GSERIALIZED_P(1);
2221 
2222  errorIfGeometryCollection(geom1,geom2);
2224 
2225  /* A.Crosses(Empty) == FALSE */
2226  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2227  PG_RETURN_BOOL(false);
2228 
2229  /*
2230  * short-circuit 1: if geom2 bounding box does not overlap
2231  * geom1 bounding box we can prematurely return FALSE.
2232  */
2233  if ( gserialized_get_gbox_p(geom1, &box1) &&
2234  gserialized_get_gbox_p(geom2, &box2) )
2235  {
2236  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2237  {
2238  PG_RETURN_BOOL(false);
2239  }
2240  }
2241 
2242  initGEOS(lwpgnotice, lwgeom_geos_error);
2243 
2244  g1 = POSTGIS2GEOS(geom1);
2245  if (!g1)
2246  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2247 
2248  g2 = POSTGIS2GEOS(geom2);
2249  if (!g2)
2250  {
2251  GEOSGeom_destroy(g1);
2253  "Second argument geometry could not be converted to GEOS");
2254  }
2255 
2256  result = GEOSCrosses(g1,g2);
2257 
2258  GEOSGeom_destroy(g1);
2259  GEOSGeom_destroy(g2);
2260 
2261  if (result == 2) HANDLE_GEOS_ERROR("GEOSCrosses");
2262 
2263  PG_FREE_IF_COPY(geom1, 0);
2264  PG_FREE_IF_COPY(geom2, 1);
2265 
2266  PG_RETURN_BOOL(result);
2267 }
2268 
2270 Datum geos_intersects(PG_FUNCTION_ARGS)
2271 {
2272  GSERIALIZED *geom1;
2273  GSERIALIZED *geom2;
2274  int result;
2275  GBOX box1, box2;
2276  PrepGeomCache *prep_cache;
2277 
2278  geom1 = PG_GETARG_GSERIALIZED_P(0);
2279  geom2 = PG_GETARG_GSERIALIZED_P(1);
2280 
2281  errorIfGeometryCollection(geom1,geom2);
2283 
2284  /* A.Intersects(Empty) == FALSE */
2285  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2286  PG_RETURN_BOOL(false);
2287 
2288  /*
2289  * short-circuit 1: if geom2 bounding box does not overlap
2290  * geom1 bounding box we can prematurely return FALSE.
2291  */
2292  if ( gserialized_get_gbox_p(geom1, &box1) &&
2293  gserialized_get_gbox_p(geom2, &box2) )
2294  {
2295  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2296  {
2297  PG_RETURN_BOOL(false);
2298  }
2299  }
2300 
2301  /*
2302  * short-circuit 2: if the geoms are a point and a polygon,
2303  * call the point_outside_polygon function.
2304  */
2305  if ((is_point(geom1) && is_poly(geom2)) || (is_poly(geom1) && is_point(geom2)))
2306  {
2307  GSERIALIZED* gpoly = is_poly(geom1) ? geom1 : geom2;
2308  GSERIALIZED* gpoint = is_point(geom1) ? geom1 : geom2;
2309  RTREE_POLY_CACHE* cache = GetRtreeCache(fcinfo, gpoly);
2310  int retval;
2311 
2312  POSTGIS_DEBUG(3, "Point in Polygon test requested...short-circuiting.");
2313  if (gserialized_get_type(gpoint) == POINTTYPE)
2314  {
2315  LWGEOM* point = lwgeom_from_gserialized(gpoint);
2316  int pip_result = pip_short_circuit(cache, lwgeom_as_lwpoint(point), gpoly);
2317  lwgeom_free(point);
2318 
2319  retval = (pip_result != -1); /* not outside */
2320  }
2321  else if (gserialized_get_type(gpoint) == MULTIPOINTTYPE)
2322  {
2324  uint32_t i;
2325 
2326  retval = LW_FALSE;
2327  for (i = 0; i < mpoint->ngeoms; i++)
2328  {
2329  int pip_result = pip_short_circuit(cache, mpoint->geoms[i], gpoly);
2330  if (pip_result != -1) /* not outside */
2331  {
2332  retval = LW_TRUE;
2333  break;
2334  }
2335  }
2336 
2337  lwmpoint_free(mpoint);
2338  }
2339  else
2340  {
2341  /* Never get here */
2342  elog(ERROR,"Type isn't point or multipoint!");
2343  PG_RETURN_NULL();
2344  }
2345 
2346  PG_FREE_IF_COPY(geom1, 0);
2347  PG_FREE_IF_COPY(geom2, 1);
2348  PG_RETURN_BOOL(retval);
2349  }
2350 
2351  initGEOS(lwpgnotice, lwgeom_geos_error);
2352  prep_cache = GetPrepGeomCache( fcinfo, geom1, geom2 );
2353 
2354  if ( prep_cache && prep_cache->prepared_geom )
2355  {
2356  if ( prep_cache->argnum == 1 )
2357  {
2358  GEOSGeometry *g = POSTGIS2GEOS(geom2);
2359  if (!g)
2360  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2361  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2362  GEOSGeom_destroy(g);
2363  }
2364  else
2365  {
2366  GEOSGeometry *g = POSTGIS2GEOS(geom1);
2367  if (!g)
2368  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
2369  result = GEOSPreparedIntersects( prep_cache->prepared_geom, g);
2370  GEOSGeom_destroy(g);
2371  }
2372  }
2373  else
2374  {
2375  GEOSGeometry *g1;
2376  GEOSGeometry *g2;
2377  g1 = POSTGIS2GEOS(geom1);
2378  if (!g1)
2380  "First argument geometry could not be converted to "
2381  "GEOS");
2382  g2 = POSTGIS2GEOS(geom2);
2383  if (!g2)
2384  {
2385  GEOSGeom_destroy(g1);
2387  "Second argument geometry could not be converted "
2388  "to GEOS");
2389  }
2390  result = GEOSIntersects( g1, g2);
2391  GEOSGeom_destroy(g1);
2392  GEOSGeom_destroy(g2);
2393  }
2394 
2395  if (result == 2) HANDLE_GEOS_ERROR("GEOSIntersects");
2396 
2397  PG_FREE_IF_COPY(geom1, 0);
2398  PG_FREE_IF_COPY(geom2, 1);
2399 
2400  PG_RETURN_BOOL(result);
2401 }
2402 
2403 
2405 Datum touches(PG_FUNCTION_ARGS)
2406 {
2407  GSERIALIZED *geom1;
2408  GSERIALIZED *geom2;
2409  GEOSGeometry *g1, *g2;
2410  bool result;
2411  GBOX box1, box2;
2412 
2413  geom1 = PG_GETARG_GSERIALIZED_P(0);
2414  geom2 = PG_GETARG_GSERIALIZED_P(1);
2415 
2416  errorIfGeometryCollection(geom1,geom2);
2418 
2419  /* A.Touches(Empty) == FALSE */
2420  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2421  PG_RETURN_BOOL(false);
2422 
2423  /*
2424  * short-circuit 1: if geom2 bounding box does not overlap
2425  * geom1 bounding box we can prematurely return FALSE.
2426  */
2427  if ( gserialized_get_gbox_p(geom1, &box1) &&
2428  gserialized_get_gbox_p(geom2, &box2) )
2429  {
2430  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2431  {
2432  PG_RETURN_BOOL(false);
2433  }
2434  }
2435 
2436  initGEOS(lwpgnotice, lwgeom_geos_error);
2437 
2438  g1 = POSTGIS2GEOS(geom1 );
2439  if (!g1)
2440  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2441 
2442  g2 = POSTGIS2GEOS(geom2 );
2443  if (!g2)
2444  {
2445  GEOSGeom_destroy(g1);
2447  "Second argument geometry could not be converted to GEOS");
2448  }
2449 
2450  result = GEOSTouches(g1,g2);
2451 
2452  GEOSGeom_destroy(g1);
2453  GEOSGeom_destroy(g2);
2454 
2455  if (result == 2) HANDLE_GEOS_ERROR("GEOSTouches");
2456 
2457  PG_FREE_IF_COPY(geom1, 0);
2458  PG_FREE_IF_COPY(geom2, 1);
2459 
2460  PG_RETURN_BOOL(result);
2461 }
2462 
2463 
2465 Datum disjoint(PG_FUNCTION_ARGS)
2466 {
2467  GSERIALIZED *geom1;
2468  GSERIALIZED *geom2;
2469  GEOSGeometry *g1, *g2;
2470  bool result;
2471  GBOX box1, box2;
2472 
2473  geom1 = PG_GETARG_GSERIALIZED_P(0);
2474  geom2 = PG_GETARG_GSERIALIZED_P(1);
2475 
2476  errorIfGeometryCollection(geom1,geom2);
2478 
2479  /* A.Disjoint(Empty) == TRUE */
2480  if ( gserialized_is_empty(geom1) || gserialized_is_empty(geom2) )
2481  PG_RETURN_BOOL(true);
2482 
2483  /*
2484  * short-circuit 1: if geom2 bounding box does not overlap
2485  * geom1 bounding box we can prematurely return TRUE.
2486  */
2487  if ( gserialized_get_gbox_p(geom1, &box1) &&
2488  gserialized_get_gbox_p(geom2, &box2) )
2489  {
2490  if ( gbox_overlaps_2d(&box1, &box2) == LW_FALSE )
2491  {
2492  PG_RETURN_BOOL(true);
2493  }
2494  }
2495 
2496  initGEOS(lwpgnotice, lwgeom_geos_error);
2497 
2498  g1 = POSTGIS2GEOS(geom1);
2499  if (!g1)
2500  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2501 
2502  g2 = POSTGIS2GEOS(geom2);
2503  if (!g2)
2504  {
2505  GEOSGeom_destroy(g1);
2507  "Second argument geometry could not be converted to GEOS");
2508  }
2509 
2510  result = GEOSDisjoint(g1,g2);
2511 
2512  GEOSGeom_destroy(g1);
2513  GEOSGeom_destroy(g2);
2514 
2515  if (result == 2) HANDLE_GEOS_ERROR("GEOSDisjoint");
2516 
2517  PG_FREE_IF_COPY(geom1, 0);
2518  PG_FREE_IF_COPY(geom2, 1);
2519 
2520  PG_RETURN_BOOL(result);
2521 }
2522 
2523 
2525 Datum relate_pattern(PG_FUNCTION_ARGS)
2526 {
2527  GSERIALIZED *geom1;
2528  GSERIALIZED *geom2;
2529  char *patt;
2530  bool result;
2531  GEOSGeometry *g1, *g2;
2532  size_t i;
2533 
2534  geom1 = PG_GETARG_GSERIALIZED_P(0);
2535  geom2 = PG_GETARG_GSERIALIZED_P(1);
2536 
2537  /* TODO handle empty */
2538 
2539  errorIfGeometryCollection(geom1,geom2);
2541 
2542  initGEOS(lwpgnotice, lwgeom_geos_error);
2543 
2544  g1 = POSTGIS2GEOS(geom1);
2545  if (!g1)
2546  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2547  g2 = POSTGIS2GEOS(geom2);
2548  if (!g2)
2549  {
2550  GEOSGeom_destroy(g1);
2552  "Second argument geometry could not be converted to GEOS");
2553  }
2554 
2555  patt = DatumGetCString(DirectFunctionCall1(textout,
2556  PointerGetDatum(PG_GETARG_DATUM(2))));
2557 
2558  /*
2559  ** Need to make sure 't' and 'f' are upper-case before handing to GEOS
2560  */
2561  for ( i = 0; i < strlen(patt); i++ )
2562  {
2563  if ( patt[i] == 't' ) patt[i] = 'T';
2564  if ( patt[i] == 'f' ) patt[i] = 'F';
2565  }
2566 
2567  result = GEOSRelatePattern(g1,g2,patt);
2568  GEOSGeom_destroy(g1);
2569  GEOSGeom_destroy(g2);
2570  pfree(patt);
2571 
2572  if (result == 2) HANDLE_GEOS_ERROR("GEOSRelatePattern");
2573 
2574  PG_FREE_IF_COPY(geom1, 0);
2575  PG_FREE_IF_COPY(geom2, 1);
2576 
2577  PG_RETURN_BOOL(result);
2578 }
2579 
2580 
2581 
2583 Datum relate_full(PG_FUNCTION_ARGS)
2584 {
2585  GSERIALIZED *geom1;
2586  GSERIALIZED *geom2;
2587  GEOSGeometry *g1, *g2;
2588  char *relate_str;
2589  text *result;
2590  int bnr = GEOSRELATE_BNR_OGC;
2591 
2592  POSTGIS_DEBUG(2, "in relate_full()");
2593 
2594  /* TODO handle empty */
2595 
2596  geom1 = PG_GETARG_GSERIALIZED_P(0);
2597  geom2 = PG_GETARG_GSERIALIZED_P(1);
2598 
2599  if ( PG_NARGS() > 2 )
2600  {
2601  bnr = PG_GETARG_INT32(2);
2602  }
2603 
2604  errorIfGeometryCollection(geom1,geom2);
2606 
2607  initGEOS(lwpgnotice, lwgeom_geos_error);
2608 
2609  g1 = POSTGIS2GEOS(geom1 );
2610  if (!g1)
2611  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2612  g2 = POSTGIS2GEOS(geom2 );
2613  if (!g2)
2614  {
2615  GEOSGeom_destroy(g1);
2617  "Second argument geometry could not be converted to GEOS");
2618  }
2619 
2620  POSTGIS_DEBUG(3, "constructed geometries ");
2621 
2622  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g1));
2623  POSTGIS_DEBUGF(3, "%s", GEOSGeomToWKT(g2));
2624 
2625  relate_str = GEOSRelateBoundaryNodeRule(g1, g2, bnr);
2626 
2627  GEOSGeom_destroy(g1);
2628  GEOSGeom_destroy(g2);
2629 
2630  if (!relate_str) HANDLE_GEOS_ERROR("GEOSRelate");
2631 
2632  result = cstring2text(relate_str);
2633  GEOSFree(relate_str);
2634 
2635  PG_FREE_IF_COPY(geom1, 0);
2636  PG_FREE_IF_COPY(geom2, 1);
2637 
2638  PG_RETURN_TEXT_P(result);
2639 }
2640 
2641 
2643 Datum ST_Equals(PG_FUNCTION_ARGS)
2644 {
2645  GSERIALIZED *geom1;
2646  GSERIALIZED *geom2;
2647  GEOSGeometry *g1, *g2;
2648  bool result;
2649  GBOX box1, box2;
2650 
2651  geom1 = PG_GETARG_GSERIALIZED_P(0);
2652  geom2 = PG_GETARG_GSERIALIZED_P(1);
2653 
2654  errorIfGeometryCollection(geom1,geom2);
2656 
2657  /* Empty == Empty */
2658  if ( gserialized_is_empty(geom1) && gserialized_is_empty(geom2) )
2659  PG_RETURN_BOOL(true);
2660 
2661  /*
2662  * short-circuit: If geom1 and geom2 do not have the same bounding box
2663  * we can return FALSE.
2664  */
2665  if ( gserialized_get_gbox_p(geom1, &box1) &&
2666  gserialized_get_gbox_p(geom2, &box2) )
2667  {
2668  if ( gbox_same_2d_float(&box1, &box2) == LW_FALSE )
2669  {
2670  PG_RETURN_BOOL(false);
2671  }
2672  }
2673 
2674  /*
2675  * short-circuit: if geom1 and geom2 are binary-equivalent, we can return
2676  * TRUE. This is much faster than doing the comparison using GEOS.
2677  */
2678  if (VARSIZE(geom1) == VARSIZE(geom2) && !memcmp(geom1, geom2, VARSIZE(geom1))) {
2679  PG_RETURN_BOOL(true);
2680  }
2681 
2682  initGEOS(lwpgnotice, lwgeom_geos_error);
2683 
2684  g1 = POSTGIS2GEOS(geom1);
2685 
2686  if (!g1)
2687  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2688 
2689  g2 = POSTGIS2GEOS(geom2);
2690 
2691  if (!g2)
2692  {
2693  GEOSGeom_destroy(g1);
2695  "Second argument geometry could not be converted to GEOS");
2696  }
2697 
2698  result = GEOSEquals(g1,g2);
2699 
2700  GEOSGeom_destroy(g1);
2701  GEOSGeom_destroy(g2);
2702 
2703  if (result == 2) HANDLE_GEOS_ERROR("GEOSEquals");
2704 
2705  PG_FREE_IF_COPY(geom1, 0);
2706  PG_FREE_IF_COPY(geom2, 1);
2707 
2708  PG_RETURN_BOOL(result);
2709 }
2710 
2712 Datum issimple(PG_FUNCTION_ARGS)
2713 {
2714  GSERIALIZED *geom;
2715  LWGEOM *lwgeom_in;
2716  int result;
2717 
2718  POSTGIS_DEBUG(2, "issimple called");
2719 
2720  geom = PG_GETARG_GSERIALIZED_P(0);
2721 
2722  if ( gserialized_is_empty(geom) )
2723  PG_RETURN_BOOL(true);
2724 
2725  lwgeom_in = lwgeom_from_gserialized(geom);
2726  result = lwgeom_is_simple(lwgeom_in);
2727  lwgeom_free(lwgeom_in) ;
2728  PG_FREE_IF_COPY(geom, 0);
2729 
2730  if (result == -1) {
2731  PG_RETURN_NULL(); /*never get here */
2732  }
2733 
2734  PG_RETURN_BOOL(result);
2735 }
2736 
2738 Datum isring(PG_FUNCTION_ARGS)
2739 {
2740  GSERIALIZED *geom;
2741  GEOSGeometry *g1;
2742  int result;
2743 
2744  geom = PG_GETARG_GSERIALIZED_P(0);
2745 
2746  /* Empty things can't close */
2747  if ( gserialized_is_empty(geom) )
2748  PG_RETURN_BOOL(false);
2749 
2750  initGEOS(lwpgnotice, lwgeom_geos_error);
2751 
2752  g1 = POSTGIS2GEOS(geom);
2753  if (!g1)
2754  HANDLE_GEOS_ERROR("First argument geometry could not be converted to GEOS");
2755 
2756  if ( GEOSGeomTypeId(g1) != GEOS_LINESTRING )
2757  {
2758  GEOSGeom_destroy(g1);
2759  elog(ERROR, "ST_IsRing() should only be called on a linear feature");
2760  }
2761 
2762  result = GEOSisRing(g1);
2763  GEOSGeom_destroy(g1);
2764 
2765  if (result == 2) HANDLE_GEOS_ERROR("GEOSisRing");
2766 
2767  PG_FREE_IF_COPY(geom, 0);
2768  PG_RETURN_BOOL(result);
2769 }
2770 
2771 
2772 
2773 
2774 
2775 GSERIALIZED *
2776 GEOS2POSTGIS(GEOSGeom geom, char want3d)
2777 {
2778  LWGEOM *lwgeom;
2779  GSERIALIZED *result;
2780 
2781  lwgeom = GEOS2LWGEOM(geom, want3d);
2782  if ( ! lwgeom )
2783  {
2784  lwpgerror("%s: GEOS2LWGEOM returned NULL", __func__);
2785  return NULL;
2786  }
2787 
2788  POSTGIS_DEBUGF(4, "%s: GEOS2LWGEOM returned a %s", __func__, lwgeom_summary(lwgeom, 0));
2789 
2790  if ( lwgeom_needs_bbox(lwgeom) == LW_TRUE )
2791  {
2792  lwgeom_add_bbox(lwgeom);
2793  }
2794 
2795  result = geometry_serialize(lwgeom);
2796  lwgeom_free(lwgeom);
2797 
2798  return result;
2799 }
2800 
2801 /*-----=POSTGIS2GEOS= */
2802 
2803 
2804 GEOSGeometry *
2806 {
2807  GEOSGeometry *ret;
2808  LWGEOM *lwgeom = lwgeom_from_gserialized(pglwgeom);
2809  if ( ! lwgeom )
2810  {
2811  lwpgerror("POSTGIS2GEOS: unable to deserialize input");
2812  return NULL;
2813  }
2814  ret = LWGEOM2GEOS(lwgeom, 0);
2815  lwgeom_free(lwgeom);
2816  if ( ! ret )
2817  {
2818  /* lwpgerror("POSTGIS2GEOS conversion failed"); */
2819  return NULL;
2820  }
2821  return ret;
2822 }
2823 
2824 uint32_t array_nelems_not_null(ArrayType* array) {
2825  ArrayIterator iterator;
2826  Datum value;
2827  bool isnull;
2828  uint32_t nelems_not_null = 0;
2829 
2830 #if POSTGIS_PGSQL_VERSION >= 95
2831  iterator = array_create_iterator(array, 0, NULL);
2832 #else
2833  iterator = array_create_iterator(array, 0);
2834 #endif
2835  while(array_iterate(iterator, &value, &isnull) )
2836  {
2837  if (!isnull)
2838  {
2839  nelems_not_null++;
2840  }
2841  }
2842  array_free_iterator(iterator);
2843 
2844  return nelems_not_null;
2845 }
2846 
2847 /* ARRAY2LWGEOM: Converts the non-null elements of a Postgres array into a LWGEOM* array */
2848 LWGEOM** ARRAY2LWGEOM(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2849 {
2850  ArrayIterator iterator;
2851  Datum value;
2852  bool isnull;
2853  bool gotsrid = false;
2854  uint32_t i = 0;
2855 
2856  LWGEOM** lw_geoms = palloc(nelems * sizeof(LWGEOM*));
2857 
2858 #if POSTGIS_PGSQL_VERSION >= 95
2859  iterator = array_create_iterator(array, 0, NULL);
2860 #else
2861  iterator = array_create_iterator(array, 0);
2862 #endif
2863 
2864  while(array_iterate(iterator, &value, &isnull))
2865  {
2866  GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2867 
2868  if (isnull)
2869  {
2870  continue;
2871  }
2872 
2873  *is3d = *is3d || gserialized_has_z(geom);
2874 
2875  lw_geoms[i] = lwgeom_from_gserialized(geom);
2876  if (!lw_geoms[i]) /* error in creation */
2877  {
2878  lwpgerror("Geometry deserializing geometry");
2879  return NULL;
2880  }
2881  if (!gotsrid)
2882  {
2883  gotsrid = true;
2884  *srid = gserialized_get_srid(geom);
2885  }
2886  else if (*srid != gserialized_get_srid(geom))
2887  {
2889  return NULL;
2890  }
2891 
2892  i++;
2893  }
2894 
2895  return lw_geoms;
2896 }
2897 
2898 /* ARRAY2GEOS: Converts the non-null elements of a Postgres array into a GEOSGeometry* array */
2899 GEOSGeometry** ARRAY2GEOS(ArrayType* array, uint32_t nelems, int* is3d, int* srid)
2900 {
2901  ArrayIterator iterator;
2902  Datum value;
2903  bool isnull;
2904  bool gotsrid = false;
2905  uint32_t i = 0;
2906 
2907  GEOSGeometry** geos_geoms = palloc(nelems * sizeof(GEOSGeometry*));
2908 
2909 #if POSTGIS_PGSQL_VERSION >= 95
2910  iterator = array_create_iterator(array, 0, NULL);
2911 #else
2912  iterator = array_create_iterator(array, 0);
2913 #endif
2914 
2915  while(array_iterate(iterator, &value, &isnull))
2916  {
2917  GSERIALIZED *geom = (GSERIALIZED*) DatumGetPointer(value);
2918 
2919  if (isnull)
2920  {
2921  continue;
2922  }
2923 
2924  *is3d = *is3d || gserialized_has_z(geom);
2925 
2926  geos_geoms[i] = POSTGIS2GEOS(geom);
2927  if (!geos_geoms[i])
2928  {
2929  uint32_t j;
2930  lwpgerror("Geometry could not be converted to GEOS");
2931 
2932  for (j = 0; j < i; j++) {
2933  GEOSGeom_destroy(geos_geoms[j]);
2934  }
2935  return NULL;
2936  }
2937 
2938  if (!gotsrid)
2939  {
2940  *srid = gserialized_get_srid(geom);
2941  gotsrid = true;
2942  }
2943  else if (*srid != gserialized_get_srid(geom))
2944  {
2945  uint32_t j;
2947 
2948  for (j = 0; j <= i; j++) {
2949  GEOSGeom_destroy(geos_geoms[j]);
2950  }
2951  return NULL;
2952  }
2953 
2954  i++;
2955  }
2956 
2957  array_free_iterator(iterator);
2958  return geos_geoms;
2959 }
2960 
2962 Datum GEOSnoop(PG_FUNCTION_ARGS)
2963 {
2964  GSERIALIZED *geom;
2965  GEOSGeometry *geosgeom;
2966  GSERIALIZED *lwgeom_result;
2967 
2968  initGEOS(lwpgnotice, lwgeom_geos_error);
2969 
2970  geom = PG_GETARG_GSERIALIZED_P(0);
2971  geosgeom = POSTGIS2GEOS(geom);
2972  if ( ! geosgeom ) PG_RETURN_NULL();
2973 
2974  lwgeom_result = GEOS2POSTGIS(geosgeom, gserialized_has_z(geom));
2975  GEOSGeom_destroy(geosgeom);
2976 
2977  PG_FREE_IF_COPY(geom, 0);
2978 
2979  PG_RETURN_POINTER(lwgeom_result);
2980 }
2981 
2983 Datum polygonize_garray(PG_FUNCTION_ARGS)
2984 {
2985  ArrayType *array;
2986  int is3d = 0;
2987  uint32 nelems, i;
2988  GSERIALIZED *result;
2989  GEOSGeometry *geos_result;
2990  const GEOSGeometry **vgeoms;
2991  int srid=SRID_UNKNOWN;
2992 #if POSTGIS_DEBUG_LEVEL >= 3
2993  static int call=1;
2994 #endif
2995 
2996 #if POSTGIS_DEBUG_LEVEL >= 3
2997  call++;
2998 #endif
2999 
3000  if (PG_ARGISNULL(0))
3001  PG_RETURN_NULL();
3002 
3003  array = PG_GETARG_ARRAYTYPE_P(0);
3004  nelems = array_nelems_not_null(array);
3005 
3006  if (nelems == 0)
3007  PG_RETURN_NULL();
3008 
3009  POSTGIS_DEBUGF(3, "polygonize_garray: number of non-null elements: %d", nelems);
3010 
3011  /* Ok, we really need geos now ;) */
3012  initGEOS(lwpgnotice, lwgeom_geos_error);
3013 
3014  vgeoms = (const GEOSGeometry**) ARRAY2GEOS(array, nelems, &is3d, &srid);
3015 
3016  POSTGIS_DEBUG(3, "polygonize_garray: invoking GEOSpolygonize");
3017 
3018  geos_result = GEOSPolygonize(vgeoms, nelems);
3019 
3020  POSTGIS_DEBUG(3, "polygonize_garray: GEOSpolygonize returned");
3021 
3022  for (i=0; i<nelems; ++i) GEOSGeom_destroy((GEOSGeometry *)vgeoms[i]);
3023  pfree(vgeoms);
3024 
3025  if ( ! geos_result ) PG_RETURN_NULL();
3026 
3027  GEOSSetSRID(geos_result, srid);
3028  result = GEOS2POSTGIS(geos_result, is3d);
3029  GEOSGeom_destroy(geos_result);
3030  if (!result)
3031  {
3032  elog(ERROR, "%s returned an error", __func__);
3033  PG_RETURN_NULL(); /*never get here */
3034  }
3035 
3036  PG_RETURN_POINTER(result);
3037 }
3038 
3039 
3041 Datum clusterintersecting_garray(PG_FUNCTION_ARGS)
3042 {
3043  Datum* result_array_data;
3044  ArrayType *array, *result;
3045  int is3d = 0;
3046  uint32 nelems, nclusters, i;
3047  GEOSGeometry **geos_inputs, **geos_results;
3048  int srid=SRID_UNKNOWN;
3049 
3050  /* Parameters used to construct a result array */
3051  int16 elmlen;
3052  bool elmbyval;
3053  char elmalign;
3054 
3055  /* Null array, null geometry (should be empty?) */
3056  if (PG_ARGISNULL(0))
3057  PG_RETURN_NULL();
3058 
3059  array = PG_GETARG_ARRAYTYPE_P(0);
3060  nelems = array_nelems_not_null(array);
3061 
3062  POSTGIS_DEBUGF(3, "clusterintersecting_garray: number of non-null elements: %d", nelems);
3063 
3064  if ( nelems == 0 ) PG_RETURN_NULL();
3065 
3066  /* TODO short-circuit for one element? */
3067 
3068  /* Ok, we really need geos now ;) */
3069  initGEOS(lwpgnotice, lwgeom_geos_error);
3070 
3071  geos_inputs = ARRAY2GEOS(array, nelems, &is3d, &srid);
3072  if(!geos_inputs)
3073  {
3074  PG_RETURN_NULL();
3075  }
3076 
3077  if (cluster_intersecting(geos_inputs, nelems, &geos_results, &nclusters) != LW_SUCCESS)
3078  {
3079  elog(ERROR, "clusterintersecting: Error performing clustering");
3080  PG_RETURN_NULL();
3081  }
3082  pfree(geos_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3083 
3084  if (!geos_results) PG_RETURN_NULL();
3085 
3086  result_array_data = palloc(nclusters * sizeof(Datum));
3087  for (i=0; i<nclusters; ++i)
3088  {
3089  result_array_data[i] = PointerGetDatum(GEOS2POSTGIS(geos_results[i], is3d));
3090  GEOSGeom_destroy(geos_results[i]);
3091  }
3092  pfree(geos_results);
3093 
3094  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3095  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3096 
3097  if (!result)
3098  {
3099  elog(ERROR, "clusterintersecting: Error constructing return-array");
3100  PG_RETURN_NULL();
3101  }
3102 
3103  PG_RETURN_POINTER(result);
3104 }
3105 
3107 Datum cluster_within_distance_garray(PG_FUNCTION_ARGS)
3108 {
3109  Datum* result_array_data;
3110  ArrayType *array, *result;
3111  int is3d = 0;
3112  uint32 nelems, nclusters, i;
3113  LWGEOM** lw_inputs;
3114  LWGEOM** lw_results;
3115  double tolerance;
3116  int srid=SRID_UNKNOWN;
3117 
3118  /* Parameters used to construct a result array */
3119  int16 elmlen;
3120  bool elmbyval;
3121  char elmalign;
3122 
3123  /* Null array, null geometry (should be empty?) */
3124  if (PG_ARGISNULL(0))
3125  PG_RETURN_NULL();
3126 
3127  array = PG_GETARG_ARRAYTYPE_P(0);
3128 
3129  tolerance = PG_GETARG_FLOAT8(1);
3130  if (tolerance < 0)
3131  {
3132  lwpgerror("Tolerance must be a positive number.");
3133  PG_RETURN_NULL();
3134  }
3135 
3136  nelems = array_nelems_not_null(array);
3137 
3138  POSTGIS_DEBUGF(3, "cluster_within_distance_garray: number of non-null elements: %d", nelems);
3139 
3140  if ( nelems == 0 ) PG_RETURN_NULL();
3141 
3142  /* TODO short-circuit for one element? */
3143 
3144  /* Ok, we really need geos now ;) */
3145  initGEOS(lwpgnotice, lwgeom_geos_error);
3146 
3147  lw_inputs = ARRAY2LWGEOM(array, nelems, &is3d, &srid);
3148  if (!lw_inputs)
3149  {
3150  PG_RETURN_NULL();
3151  }
3152 
3153  if (cluster_within_distance(lw_inputs, nelems, tolerance, &lw_results, &nclusters) != LW_SUCCESS)
3154  {
3155  elog(ERROR, "cluster_within: Error performing clustering");
3156  PG_RETURN_NULL();
3157  }
3158  pfree(lw_inputs); /* don't need to destroy items because GeometryCollections have taken ownership */
3159 
3160  if (!lw_results) PG_RETURN_NULL();
3161 
3162  result_array_data = palloc(nclusters * sizeof(Datum));
3163  for (i=0; i<nclusters; ++i)
3164  {
3165  result_array_data[i] = PointerGetDatum(gserialized_from_lwgeom(lw_results[i], NULL));
3166  lwgeom_free(lw_results[i]);
3167  }
3168  pfree(lw_results);
3169 
3170  get_typlenbyvalalign(array->elemtype, &elmlen, &elmbyval, &elmalign);
3171  result = construct_array(result_array_data, nclusters, array->elemtype, elmlen, elmbyval, elmalign);
3172 
3173  if (!result)
3174  {
3175  elog(ERROR, "clusterwithin: Error constructing return-array");
3176  PG_RETURN_NULL();
3177  }
3178 
3179  PG_RETURN_POINTER(result);
3180 }
3181 
3183 Datum linemerge(PG_FUNCTION_ARGS)
3184 {
3185  GSERIALIZED *geom1;
3186  GSERIALIZED *result;
3187  LWGEOM *lwgeom1, *lwresult ;
3188 
3189  geom1 = PG_GETARG_GSERIALIZED_P(0);
3190 
3191 
3192  lwgeom1 = lwgeom_from_gserialized(geom1) ;
3193 
3194  lwresult = lwgeom_linemerge(lwgeom1);
3195  result = geometry_serialize(lwresult) ;
3196 
3197  lwgeom_free(lwgeom1) ;
3198  lwgeom_free(lwresult) ;
3199 
3200  PG_FREE_IF_COPY(geom1, 0);
3201 
3202  PG_RETURN_POINTER(result);
3203 }
3204 
3205 /*
3206  * Take a geometry and return an areal geometry
3207  * (Polygon or MultiPolygon).
3208  * Actually a wrapper around GEOSpolygonize,
3209  * transforming the resulting collection into
3210  * a valid polygon Geometry.
3211  */
3213 Datum ST_BuildArea(PG_FUNCTION_ARGS)
3214 {
3215  GSERIALIZED *result;
3216  GSERIALIZED *geom;
3217  LWGEOM *lwgeom_in, *lwgeom_out;
3218 
3219  geom = PG_GETARG_GSERIALIZED_P(0);
3220  lwgeom_in = lwgeom_from_gserialized(geom);
3221 
3222  lwgeom_out = lwgeom_buildarea(lwgeom_in);
3223  lwgeom_free(lwgeom_in) ;
3224 
3225  if ( ! lwgeom_out )
3226  {
3227  PG_FREE_IF_COPY(geom, 0);
3228  PG_RETURN_NULL();
3229  }
3230 
3231  result = geometry_serialize(lwgeom_out) ;
3232  lwgeom_free(lwgeom_out) ;
3233 
3234  PG_FREE_IF_COPY(geom, 0);
3235  PG_RETURN_POINTER(result);
3236 }
3237 
3238 /*
3239  * Take the vertices of a geometry and builds
3240  * Delaunay triangles around them.
3241  */
3243 Datum ST_DelaunayTriangles(PG_FUNCTION_ARGS)
3244 {
3245  GSERIALIZED *result;
3246  GSERIALIZED *geom;
3247  LWGEOM *lwgeom_in, *lwgeom_out;
3248  double tolerance = 0.0;
3249  int flags = 0;
3250 
3251  geom = PG_GETARG_GSERIALIZED_P(0);
3252  tolerance = PG_GETARG_FLOAT8(1);
3253  flags = PG_GETARG_INT32(2);
3254 
3255  lwgeom_in = lwgeom_from_gserialized(geom);
3256  lwgeom_out = lwgeom_delaunay_triangulation(lwgeom_in, tolerance, flags);
3257  lwgeom_free(lwgeom_in) ;
3258 
3259  if ( ! lwgeom_out )
3260  {
3261  PG_FREE_IF_COPY(geom, 0);
3262  PG_RETURN_NULL();
3263  }
3264 
3265  result = geometry_serialize(lwgeom_out) ;
3266  lwgeom_free(lwgeom_out) ;
3267 
3268  PG_FREE_IF_COPY(geom, 0);
3269  PG_RETURN_POINTER(result);
3270 }
3271 
3272 /*
3273  * ST_Snap
3274  *
3275  * Snap a geometry to another with a given tolerance
3276  */
3277 Datum ST_Snap(PG_FUNCTION_ARGS);
3279 Datum ST_Snap(PG_FUNCTION_ARGS)
3280 {
3281  GSERIALIZED *geom1, *geom2, *result;
3282  LWGEOM *lwgeom1, *lwgeom2, *lwresult;
3283  double tolerance;
3284 
3285  geom1 = PG_GETARG_GSERIALIZED_P(0);
3286  geom2 = PG_GETARG_GSERIALIZED_P(1);
3287  tolerance = PG_GETARG_FLOAT8(2);
3288 
3289  lwgeom1 = lwgeom_from_gserialized(geom1);
3290  lwgeom2 = lwgeom_from_gserialized(geom2);
3291 
3292  lwresult = lwgeom_snap(lwgeom1, lwgeom2, tolerance);
3293  lwgeom_free(lwgeom1);
3294  lwgeom_free(lwgeom2);
3295  PG_FREE_IF_COPY(geom1, 0);
3296  PG_FREE_IF_COPY(geom2, 1);
3297 
3298  result = geometry_serialize(lwresult);
3299  lwgeom_free(lwresult);
3300 
3301  PG_RETURN_POINTER(result);
3302 }
3303 
3304 /*
3305  * ST_Split
3306  *
3307  * Split polygon by line, line by line, line by point.
3308  * Returns at most components as a collection.
3309  * First element of the collection is always the part which
3310  * remains after the cut, while the second element is the
3311  * part which has been cut out. We arbitrarely take the part
3312  * on the *right* of cut lines as the part which has been cut out.
3313  * For a line cut by a point the part which remains is the one
3314  * from start of the line to the cut point.
3315  *
3316  *
3317  * Author: Sandro Santilli <strk@kbt.io>
3318  *
3319  * Work done for Faunalia (http://www.faunalia.it) with fundings
3320  * from Regione Toscana - Sistema Informativo per il Governo
3321  * del Territorio e dell'Ambiente (RT-SIGTA).
3322  *
3323  * Thanks to the PostGIS community for sharing poly/line ideas [1]
3324  *
3325  * [1] http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithLineString
3326  *
3327  */
3328 Datum ST_Split(PG_FUNCTION_ARGS);
3330 Datum ST_Split(PG_FUNCTION_ARGS)
3331 {
3332  GSERIALIZED *in, *blade_in, *out;
3333  LWGEOM *lwgeom_in, *lwblade_in, *lwgeom_out;
3334 
3335  in = PG_GETARG_GSERIALIZED_P(0);
3336  lwgeom_in = lwgeom_from_gserialized(in);
3337 
3338  blade_in = PG_GETARG_GSERIALIZED_P(1);
3339  lwblade_in = lwgeom_from_gserialized(blade_in);
3340 
3341  error_if_srid_mismatch(lwgeom_in->srid, lwblade_in->srid);
3342 
3343  lwgeom_out = lwgeom_split(lwgeom_in, lwblade_in);
3344  lwgeom_free(lwgeom_in);
3345  lwgeom_free(lwblade_in);
3346 
3347  if ( ! lwgeom_out )
3348  {
3349  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3350  PG_FREE_IF_COPY(blade_in, 1);
3351  PG_RETURN_NULL();
3352  }
3353 
3354  out = geometry_serialize(lwgeom_out);
3355  lwgeom_free(lwgeom_out);
3356  PG_FREE_IF_COPY(in, 0); /* possibly referenced by lwgeom_out */
3357  PG_FREE_IF_COPY(blade_in, 1);
3358 
3359  PG_RETURN_POINTER(out);
3360 }
3361 
3362 /**********************************************************************
3363  *
3364  * ST_SharedPaths
3365  *
3366  * Return the set of paths shared between two linear geometries,
3367  * and their direction (same or opposite).
3368  *
3369  * Developed by Sandro Santilli (strk@kbt.io) for Faunalia
3370  * (http://www.faunalia.it) with funding from Regione Toscana - Sistema
3371  * Informativo per la Gestione del Territorio e dell' Ambiente
3372  * [RT-SIGTA]". For the project: "Sviluppo strumenti software per il
3373  * trattamento di dati geografici basati su QuantumGIS e Postgis (CIG
3374  * 0494241492)"
3375  *
3376  **********************************************************************/
3377 Datum ST_SharedPaths(PG_FUNCTION_ARGS);
3379 Datum ST_SharedPaths(PG_FUNCTION_ARGS)
3380 {
3381  GSERIALIZED *geom1, *geom2, *out;
3382  LWGEOM *g1, *g2, *lwgeom_out;
3383 
3384  geom1 = PG_GETARG_GSERIALIZED_P(0);
3385  geom2 = PG_GETARG_GSERIALIZED_P(1);
3386 
3387  g1 = lwgeom_from_gserialized(geom1);
3388  g2 = lwgeom_from_gserialized(geom2);
3389 
3390  lwgeom_out = lwgeom_sharedpaths(g1, g2);
3391  lwgeom_free(g1);
3392  lwgeom_free(g2);
3393 
3394  if ( ! lwgeom_out )
3395  {
3396  PG_FREE_IF_COPY(geom1, 0);
3397  PG_FREE_IF_COPY(geom2, 1);
3398  PG_RETURN_NULL();
3399  }
3400 
3401  out = geometry_serialize(lwgeom_out);
3402  lwgeom_free(lwgeom_out);
3403 
3404  PG_FREE_IF_COPY(geom1, 0);
3405  PG_FREE_IF_COPY(geom2, 1);
3406  PG_RETURN_POINTER(out);
3407 }
3408 
3409 
3410 /**********************************************************************
3411  *
3412  * ST_Node
3413  *
3414  * Fully node a set of lines using the least possible nodes while
3415  * preserving all of the input ones.
3416  *
3417  **********************************************************************/
3418 Datum ST_Node(PG_FUNCTION_ARGS);
3420 Datum ST_Node(PG_FUNCTION_ARGS)
3421 {
3422  GSERIALIZED *geom1, *out;
3423  LWGEOM *g1, *lwgeom_out;
3424 
3425  geom1 = PG_GETARG_GSERIALIZED_P(0);
3426 
3427  g1 = lwgeom_from_gserialized(geom1);
3428 
3429  lwgeom_out = lwgeom_node(g1);
3430  lwgeom_free(g1);
3431 
3432  if ( ! lwgeom_out )
3433  {
3434  PG_FREE_IF_COPY(geom1, 0);
3435  PG_RETURN_NULL();
3436  }
3437 
3438  out = geometry_serialize(lwgeom_out);
3439  lwgeom_free(lwgeom_out);
3440 
3441  PG_FREE_IF_COPY(geom1, 0);
3442  PG_RETURN_POINTER(out);
3443 }
3444 
3445 /******************************************
3446  *
3447  * ST_Voronoi
3448  *
3449  * Returns a Voronoi diagram constructed
3450  * from the points of the input geometry.
3451  *
3452  ******************************************/
3453 Datum ST_Voronoi(PG_FUNCTION_ARGS);
3455 Datum ST_Voronoi(PG_FUNCTION_ARGS)
3456 {
3457 #if POSTGIS_GEOS_VERSION < 35
3458  lwpgerror("The GEOS version this PostGIS binary "
3459  "was compiled against (%d) doesn't support "
3460  "'ST_Voronoi' function (3.5.0+ required)",
3462  PG_RETURN_NULL();
3463 #else /* POSTGIS_GEOS_VERSION >= 35 */
3464  GSERIALIZED* input;
3465  GSERIALIZED* clip;
3466  GSERIALIZED* result;
3467  LWGEOM* lwgeom_input;
3468  LWGEOM* lwgeom_result;
3469  double tolerance;
3470  GBOX clip_envelope;
3471  int custom_clip_envelope;
3472  int return_polygons;
3473 
3474  /* Return NULL on NULL geometry */
3475  if (PG_ARGISNULL(0))
3476  PG_RETURN_NULL();
3477 
3478  /* Read our tolerance value */
3479  if (PG_ARGISNULL(2))
3480  {
3481  lwpgerror("Tolerance must be a positive number.");
3482  PG_RETURN_NULL();
3483  }
3484 
3485  tolerance = PG_GETARG_FLOAT8(2);
3486 
3487  if (tolerance < 0)
3488  {
3489  lwpgerror("Tolerance must be a positive number.");
3490  PG_RETURN_NULL();
3491  }
3492 
3493  /* Are we returning lines or polygons? */
3494  if (PG_ARGISNULL(3))
3495  {
3496  lwpgerror("return_polygons must be true or false.");
3497  PG_RETURN_NULL();
3498  }
3499  return_polygons = PG_GETARG_BOOL(3);
3500 
3501  /* Read our clipping envelope, if applicable. */
3502  custom_clip_envelope = !PG_ARGISNULL(1);
3503  if (custom_clip_envelope) {
3504  clip = PG_GETARG_GSERIALIZED_P(1);
3505  if (!gserialized_get_gbox_p(clip, &clip_envelope))
3506  {
3507  lwpgerror("Could not determine envelope of clipping geometry.");
3508  PG_FREE_IF_COPY(clip, 1);
3509  PG_RETURN_NULL();
3510  }
3511  PG_FREE_IF_COPY(clip, 1);
3512  }
3513 
3514  /* Read our input geometry */
3515  input = PG_GETARG_GSERIALIZED_P(0);
3516 
3517  lwgeom_input = lwgeom_from_gserialized(input);
3518 
3519  if(!lwgeom_input)
3520  {
3521  lwpgerror("Could not read input geometry.");
3522  PG_FREE_IF_COPY(input, 0);
3523  PG_RETURN_NULL();
3524  }
3525 
3526  lwgeom_result = lwgeom_voronoi_diagram(lwgeom_input, custom_clip_envelope ? &clip_envelope : NULL, tolerance, !return_polygons);
3527  lwgeom_free(lwgeom_input);
3528 
3529  if (!lwgeom_result)
3530  {
3531  lwpgerror("Error computing Voronoi diagram.");
3532  PG_FREE_IF_COPY(input, 0);
3533  PG_RETURN_NULL();
3534  }
3535 
3536  result = geometry_serialize(lwgeom_result);
3537  lwgeom_free(lwgeom_result);
3538 
3539  PG_FREE_IF_COPY(input, 0);
3540  PG_RETURN_POINTER(result);
3541 
3542 #endif /* POSTGIS_GEOS_VERSION >= 35 */
3543 }
3544 
3545 /******************************************
3546  *
3547  * ST_MinimumClearance
3548  *
3549  * Returns the minimum clearance of a geometry.
3550  *
3551  ******************************************/
3552 Datum ST_MinimumClearance(PG_FUNCTION_ARGS);
3554 Datum ST_MinimumClearance(PG_FUNCTION_ARGS)
3555 {
3556 #if POSTGIS_GEOS_VERSION < 36
3557  lwpgerror("The GEOS version this PostGIS binary "
3558  "was compiled against (%d) doesn't support "
3559  "'ST_MinimumClearance' function (3.6.0+ required)",
3561  PG_RETURN_NULL();
3562 #else /* POSTGIS_GEOS_VERSION >= 36 */
3563  GSERIALIZED* input;
3564  GEOSGeometry* input_geos;
3565  int error;
3566  double result;
3567 
3568  initGEOS(lwpgnotice, lwgeom_geos_error);
3569 
3570  input = PG_GETARG_GSERIALIZED_P(0);
3571  input_geos = POSTGIS2GEOS(input);
3572  if (!input_geos)
3573  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3574 
3575  error = GEOSMinimumClearance(input_geos, &result);
3576  GEOSGeom_destroy(input_geos);
3577  if (error) HANDLE_GEOS_ERROR("Error computing minimum clearance");
3578 
3579  PG_FREE_IF_COPY(input, 0);
3580  PG_RETURN_FLOAT8(result);
3581 #endif
3582 }
3583 
3584 /******************************************
3585  *
3586  * ST_MinimumClearanceLine
3587  *
3588  * Returns the minimum clearance line of a geometry.
3589  *
3590  ******************************************/
3591 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS);
3593 Datum ST_MinimumClearanceLine(PG_FUNCTION_ARGS)
3594 {
3595 #if POSTGIS_GEOS_VERSION < 36
3596  lwpgerror("The GEOS version this PostGIS binary "
3597  "was compiled against (%d) doesn't support "
3598  "'ST_MinimumClearanceLine' function (3.6.0+ required)",
3600  PG_RETURN_NULL();
3601 #else /* POSTGIS_GEOS_VERSION >= 36 */
3602  GSERIALIZED* input;
3603  GSERIALIZED* result;
3604  GEOSGeometry* input_geos;
3605  GEOSGeometry* result_geos;
3606  int srid;
3607 
3608  initGEOS(lwpgnotice, lwgeom_geos_error);
3609 
3610  input = PG_GETARG_GSERIALIZED_P(0);
3611  srid = gserialized_get_srid(input);
3612  input_geos = POSTGIS2GEOS(input);
3613  if (!input_geos)
3614  HANDLE_GEOS_ERROR("Geometry could not be converted to GEOS");
3615 
3616  result_geos = GEOSMinimumClearanceLine(input_geos);
3617  GEOSGeom_destroy(input_geos);
3618  if (!result_geos)
3619  HANDLE_GEOS_ERROR("Error computing minimum clearance");
3620 
3621  GEOSSetSRID(result_geos, srid);
3622  result = GEOS2POSTGIS(result_geos, LW_FALSE);
3623  GEOSGeom_destroy(result_geos);
3624 
3625  PG_FREE_IF_COPY(input, 0);
3626  PG_RETURN_POINTER(result);
3627 #endif
3628 }
3629 
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
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:714
unsigned int int32
Definition: shpopen.c:273
#define POSTGIS_GEOS_VERSION
Definition: sqldefines.h:10
GBOX * bbox
Definition: liblwgeom.h:397
uint32_t ngeoms
Definition: liblwgeom.h:467
Datum covers(PG_FUNCTION_ARGS)
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:916
Datum ST_SharedPaths(PG_FUNCTION_ARGS)
#define POLYGONTYPE
Definition: liblwgeom.h:86
uint8_t flags
Definition: liblwgeom.h:396
Datum ST_Equals(PG_FUNCTION_ARGS)
Datum postgis_geos_version(PG_FUNCTION_ARGS)
double xmax
Definition: liblwgeom.h:292
Datum ST_UnaryUnion(PG_FUNCTION_ARGS)
void lwpoint_free(LWPOINT *pt)
Definition: lwpoint.c:213
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1137
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
The tree structure used for fast P-i-P tests by point_in_multipolygon_rtree()
Definition: lwgeom_rtree.h:57
Datum isvalid(PG_FUNCTION_ARGS)
LWGEOM * lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyle, double mitreLimit)
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)
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:205
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:160
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:398
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:319
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:250
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:293
PG_FUNCTION_INFO_V1(postgis_geos_version)
LWGEOM * geom
double xmin
Definition: liblwgeom.h:291
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
LWGEOM * lwgeom_delaunay_triangulation(const LWGEOM *geom, double tolerance, int edgeOnly)
Take vertices of a geometry and build a delaunay triangulation on them.
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:187
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
char * text2cstring(const text *textptr)
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:727
int count
Definition: genraster.py:56
LWPOINT ** geoms
Definition: liblwgeom.h:469
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:294
Datum convexhull(PG_FUNCTION_ARGS)
Datum containsproperly(PG_FUNCTION_ARGS)
GEOSGeometry * LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix)
LWGEOM * lwgeom_construct_empty(uint8_t type, int srid, char hasz, char hasm)
Definition: lwgeom.c:2078
#define WKB_EXTENDED
Definition: liblwgeom.h:2061
Datum disjoint(PG_FUNCTION_ARGS)
Datum ST_ClipByBox2d(PG_FUNCTION_ARGS)
LWLINE * lwgeom_as_lwline(const LWGEOM *lwgeom)
Definition: lwgeom.c:169
uint8_t flags
Definition: liblwgeom.h:290
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:2069
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:159
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:232
void lwgeom_add_bbox(LWGEOM *lwgeom)
Compute a bbox if not already computed.
Definition: lwgeom.c:679
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
Datum ST_GeneratePoints(PG_FUNCTION_ARGS)
uint8_t type
Definition: liblwgeom.h:395
PrepGeomCache * GetPrepGeomCache(FunctionCallInfoData *fcinfo, GSERIALIZED *g1, GSERIALIZED *g2)
Given a couple potential geometries and a function call context, return a prepared structure for one ...
Datum geos_difference(PG_FUNCTION_ARGS)
GSERIALIZED * GEOS2POSTGIS(GEOSGeom geom, char want3d)
LWGEOM * lwpoint_as_lwgeom(const LWPOINT *obj)
Definition: lwgeom.c:334
Datum isvaliddetail(PG_FUNCTION_ARGS)
#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:59
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:1386
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:1190
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.