source: trunk/postgis/lwgeom_geos.c

Last change on this file was 17951, checked in by Raul Marin, 5 years ago

Use get_call_result_type to retrieve tuple descriptions

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