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